MATLAB: Conditional statements and ODEs

functionodetracer study

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?
%

Best Answer

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.
Related Question