MATLAB: Looping ODEs, Attempted to access Cu(2); index out of bounds because numel(Cu)=1.

accessattemptedindexloopnumelodeto

clc;
Funderin(1) = 1.7153;
Funderout(1) = 1;
Fsep(1) = 0.9406;
MassUnders(1) = 12.0069;
tspan1 = linspace(0, 3, 25);
tspan2 = linspace(3, 5, 25);
Funderin(2) = 1;
Funderout(2) = 0.5;
Fsep(2) = 0.2;
MassUnders(2) = 8;
Funderin(3) = 0.5;
Funderout(3) = 0.2;
Fsep(3) = 0.1;
MassUnders(3) = 5;
Cu_in(1) = 0.9717   %tracer flow
Cui = 0;        %initial concentration
rhs = @(t,Cu,Cu_in) (Funderin(1)* Cu_in(1) - Funderout(1)*Cu(1) - Fsep(1)*Cu(1))/MassUnders(1);
rhs = @(t,Cu,Cu_in) (Funderin(2)* Cu(1) - Funderout(2)*Cu(2) - Fsep(2)*Cu(2))/MassUnders(2);
rhs = @(t,Cu,Cu_in) (Funderin(3)* Cu(2) - Funderout(3)*Cu(3) - Fsep(3)*Cu(3))/MassUnders(3);
[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('Time'); ylabel('Cu');
The end goal is to get rhs looped for i = 1:K, where K can be upwards 50. But I can't even seem to get Cu to 'save' properly. Any ideas?

Best Answer

Now that I understand how you want the loop to work, I can help you with the code. I changed ‘rhs’ to accommodate all the parameters you want to vary, and I created vectors for them. It doesn’t produce the same plot yesterday’s code did, but then the parameters are different. The code runs without error. To run more segments, increase the parameter vector lengths. They must all be the same lengths. The rest of the code (the loop) will adapt to the parameter vector lengths, so no other changes are necessary in this code.
Change the second argument in the ‘tv’ vector linspace call if you want time spans of other than 3 time units (a guess on my part, from yesterday’s Question). The length of 25 seems to give an acceptably smooth curve, and since you may want to iterate this 50 times, gives a manageable row size for ‘cu_int’ of 1250.
See if this does what you want:
Funderin = [1.7153, 1, 0.5];
Funderout = [1, 0.5, 0.2];
Fsep = [0.9406, 0.2, 0.1];
MassUnders = [12.0069, 8, 5];
tv = linspace(0, 3, 25);
t = 0;
Cu_in = 0.9717   %tracer flow
Cui = Cu_in;        %initial concentration
Cu_int = [];
rhs = @(t,Cu,Cu_in,Fin,Fout,Fsp,MassU) (Fin* Cu_in - Fout*Cu - Fsp*Cu)/MassU;
for k1 = 1:length(Funderin)
    tspan = tv + t(end);
    [t,Cu] = ode45( @(t,Cu) rhs(t,Cu,Cu_in,Funderin(k1),Funderout(k1),Fsep(k1),MassUnders(k1)), tspan, Cui);
    Cu_in = Cu(end);
    Cui = Cu(end);
    Cu_int = [Cu_int; t, Cu];
end
plot(Cu_int(:,1),Cu_int(:,2))
xlabel('Time (Units)'); ylabel('Cu (Units)')
Related Question