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.
