Conditional statements and ODEs

function first_oder_ode

Funderin = 1.7153
Funderout = 1.7153
Fsep= 0.9406
MassUnders = 12.0069

t=0:0.01:5;

if 0 < t < 3
        Cu_in = 0.9717
else
        Cu_in = 0 
end

      [t,Cu]=ode45( @rhs, t, 0);
      plot(t,Cu);
      xlabel('t'); ylabel('Cu');

      function dCudt = rhs(t,Cu)
      dCudt = (Funderin* Cu_in - Funderout*Cu - Fsep*Cu)/MassUnders;
      end
  end

Currently with the code above I'm trying to generate a graph that should look like this with the red line instead of the blue line. It's a conditional statement that turns off Cu_in (after 3 seconds, Cu_in = 0) and what's left of it is 'drained' away by Fsep and Funderout. Have I approached the question correctly or am I missing something here?

%

Answer 1

[Star Strider]

Numerical ODE solvers don’t like discontinuities, so you have to ‘break’ the integration into two parts, using the last output of the first integration as the initial condition for the second. I coded ‘rhs’ as an anonymous function because it’s easier.

See if this does what you want:

Funderin = 1.7153;
Funderout = 1.7153;
Fsep= 0.9406;
MassUnders = 12.0069;

rhs = @(t,Cu,Cu_in) (Funderin* Cu_in - Funderout*Cu - Fsep*Cu)/MassUnders;

N = 25;
tspan1 = linspace(0, 3, N);
tspan2 = linspace(3, 5, N);

Cu_in = 0.9717
Cui = 0;
[t,Cu]=ode45( @(t,Cu) rhs(t,Cu,Cu_in), tspan1, Cui);

Cu_int = [t, Cu];

Cu_in = 0;
Cui = Cu(end);
[t,Cu]=ode45( @(t,Cu) rhs(t,Cu,Cu_in), tspan2, Cui);

Cu_int = [Cu_int; t, Cu];

plot(Cu_int(:,1),Cu_int(:,2));
xlabel('t'); ylabel('Cu');

EDIT — Forgot to include the plot call. Added it.

[SuppoCheng]

Perfect. Thank you so much. I can see what you mean by MATLAB not liking discontinuities in ODEs!

[Star Strider]

My pleasure!