for loop within an odes15 function, and plot on one graph

K=30;
lambda1=1;
lambda2=2;
T_final = 10;

f1 = figure;
figure(f1);
hold on

[t,M] = ode15s(@(t,M) RHS_focus(M,K,lambda2,lambda1), [0 T_final], [0 100]);

plot(t,M(:,1:K)','k--')
function dMdt=RHS_focus(M,K,lambda1,lambda2)
dMdt = zeros(K,1);
dMdt(1) = lambda2*M(2) - lambda1*M(1);     % Note that this is for state 0.
for i= 2:K-1
    dMdt(i) = lambda1*M(i-1) - (lambda1+lambda2)*M(i) + lambda2*M(i+1);  % For state 1 to 28
end
dMdt(K) = lambda1*M(K-1) - lambda2*M(K);    % For state 29
dMdt= M';
end

ERROR: Index exceeds the number of array elements. Index must not exceed 2.

It didnt work when i didnt use a for loop either.

I want to plot 30 ODES on one graph.

Answer 1

[Star Strider]

In order for ode15s not to complain, it will be necessary to supply a vector of 30 initial conditions to the ode15s call. Even then, things might not go completely smoothly.

K=30;
lambda1=1;
lambda2=2;
T_final = 10;

ic = linspace(0, 100, K);                                                  % Initial Conditions Vector
[t,M] = ode15s(@(t,M) RHS_focus(M,K,lambda2,lambda1), [0 T_final], ic);

f1 = figure;
figure(f1);
hold on
plot(t,M(:,1:K)','k--')
hold off

function dMdt=RHS_focus(M,K,lambda1,lambda2)
    dMdt = zeros(K,1);
    dMdt(1,:) = lambda2*M(2) - lambda1*M(1);     % Note that this is for state 0.
    for i= 2:K-1
        dMdt(i,:) = lambda1*M(i-1) - (lambda1+lambda2)*M(i) + lambda2*M(i+1);  % For state 1 to 28
    end
    dMdt(K,:) = lambda1*M(K-1) - lambda2*M(K);    % For state 29
%     dMdt= M';
end

EDIT — (21 Dec 2022 at 17:00)

Changed ‘ic’ so that it scales with ‘K’.

.