Putting multiple ODE equations into one script

I am trying to put multiple ODEs into one script. They model a CSTR in series, where each CSTR is it's own ODE. The first equation (dvdt) is supposed to model the water flow in the module and the second equation (dbdt) models the solute flow. There are 44 CSTRs in the module (ps.N)I , so I used a for loop.

function dkdt = ode_CSTR_series(t,x, ps)

dkdt = zeros(2,ps.N);

v = x(1);
b = x(2);

dvdt(1) = (ps.v_in - v(1)) - (ps.k .* ps.area .* ps.delta_P);
dbdt(1) = (ps.b_in - b(1))/ps.tau;

for i = 2:ps.N
    dvdt(i) = (v(i-1) - v(i)) - (ps.k .* ps.area .* ps.delta_P);
    dbdt(i) = (b(i-1) - b(i))/ps.tau;
end 

dkdt = [dvdt; dbdt];
end

I am trying to get two solutions using ode15s. This is my separate script.

%Initial conditions and time span
initial = zeros(2,ps.N);
initial(1,1) = 0; % mol/s
initial(2,1) = 0;
tspan = t_water0;

[t,b] = ode15s(@ode_CSTR_series, tspan, initial, [], ps);

I get an error that says "Index exceeds the number of array elements (1)". As I understand, each equation should be kept in one row.

How do I solve this? Thanks.

Answer 1

[Alan Stevens]

In

v = x(1);
b = x(2);
dvdt(1) = (ps.v_in - v(1)) - (ps.k .* ps.area .* ps.delta_P);
dbdt(1) = (ps.b_in - b(1))/ps.tau;

v and b are just scalars, so the latter two equations shoud be

dvdt(1) = (ps.v_in - v) - (ps.k .* ps.area .* ps.delta_P);
dbdt(1) = (ps.b_in - b)/ps.tau;

This explains why you get the error message, but doesn't completely solve your problem as you need v and b to be vectors, but they must be updated before being used in your for loop. You probably need to put the loop around the calling function ([t,b] = ode15s(...).

[Zafina]

How do I make a loop around my calling function?

[Alan Stevens]

Something along the following lines perhaps:

%Initial conditions and time span
tspan = t_water0;
vb0 = [0, 0];

for k = 1:ps.N
    
    [t,vb] = ode15s(@ode_CSTR_series, tspan, vb0,[],ps);
    vb0 = vb(end,end);  % i.e. select last value of v and b
                        % as initial value for next iteration 
end


function dvbdt = ode_CSTR_series(~,vb,ps)

        v = vb(1);
        b = vb(2);

        dvdt = (ps.v_in - v) - (ps.k .* ps.area .* ps.delta_P);
        dbdt = (ps.b_in - b)/ps.tau;

        dvbdt = [dvdt; dbdt];
end

[Zafina]

In the solution, k is not within any equation of the for loop, so how does it know to keep running until ps.N?

[Alan Stevens]

When k is 1 all the calculations within the loop are carried out. Then k increments to 2 and all the calculations within the loop are carried out again. This repeats until k exceeds ps.N. This happens even if k isn't used explicitly within the loop. Incidentally, I just listed a skeleton code. I assume you will want to do other things within the loop (e.g. saving or plotting values of v an b).