I am looking to repeat my Euler's method 'for' loop for different initial conditions to produce a graph with multiple plots

clear
 
%Runge-Kutta analysis of predator/prey. 
A=1;
B=0.1;
C=1;
D=0.1;
 
t0=0; tn=10; %initial and final times
 
t(1)=t0
y(1,1)=25; y(2,1)=25; %initial conditions
 
n=200; %number of steps
h=(tn-t0)/n; %step size
Exactsol1(:,1)=[0;0];
Exactsol2(:,1)=[0;0];
 
M=[0 (-C*B)/D; (A*D)/B 0] %Jacobian Matrix
 
[v lambda]=eig(M)
 
for j=1:n %loop for RK4's formula
    t(j+1)=t0+j*h
    
    k1=h*predatorRK4(t(j),y(:,j),A,B,C,D);
    k2=h*predatorRK4(t(j)+(h/2),y(:,j)+(k1/2),A,B,C,D);
    k3=h*predatorRK4(t(j)+(h/2),y(:,j)+(k2/2),A,B,C,D);
    k4=h*predatorRK4(t(j)+h,y(:,j)+k3,A,B,C,D);
    
    Exactsol1(:,j+1) = v(:,1)*exp(lambda(1,1)*t(j+1));
    Exactsol2(:,j+1) = v(:,2)*exp(lambda(2,2)*t(j+1));
    y(:,j+1)=y(:,j)+(1/6)*(k1+2*k2+2*k3+k4);
end
figure
plot(y(1,:),y(2,:))

[Muhammad Sinan]

@Matthew Stoner You can use the following code, for every ODE just change the position of "int". While for the initial conditions, you can use any format with a for loop; like

for int = [5, 10, 15, 20, 25, 30],
    ......
end

or

for int = 5:5:30
    ......
end

or

inital_conditions = [5, 10, 15, 20, 25, 30];
for i = 1:length(inital_conditions)
    ......
end

While, I have just used the first case in the following code.

clear all; close all; clf; clc;
 
%Runge-Kutta analysis of predator/prey. 
A=1;
B=0.1;
C=1;
D=0.1;
 
t0=0; tn=10; %initial and final times
 
t(1)=t0;
for int = [5, 10, 15, 20, 25, 30]
y(1,1)=int; 
y(2,1)=25; %initial conditions
 
n=200; %number of steps
h=(tn-t0)/n; %step size
Exactsol1(:,1)=[0;0];
Exactsol2(:,1)=[0;0];
 
M=[0 (-C*B)/D; (A*D)/B 0]; %Jacobian Matrix
 
[v lambda]=eig(M);
 
for j=1:n %loop for RK4's formula
    t(j+1)=t0+j*h;
    
    k1=h*predatorRK4(t(j),y(:,j),A,B,C,D);
    k2=h*predatorRK4(t(j)+(h/2),y(:,j)+(k1/2),A,B,C,D);
    k3=h*predatorRK4(t(j)+(h/2),y(:,j)+(k2/2),A,B,C,D);
    k4=h*predatorRK4(t(j)+h,y(:,j)+k3,A,B,C,D);
    
    Exactsol1(:,j+1) = v(:,1)*exp(lambda(1,1)*t(j+1));
    Exactsol2(:,j+1) = v(:,2)*exp(lambda(2,2)*t(j+1));
    y(:,j+1)=y(:,j)+(1/6)*(k1+2*k2+2*k3+k4);
end
figure
plot(y(1,:),y(2,:))
end
end

Hope, it helps you!

Answer 1

[Cris LaPierre]

You probably need to use a second for-loop then. I'll leave it to you to implement, but here's a simple example.

n = [10 200 3000]

for r = 1:length(n)
    for c = 1:n(r)
        ...
        y(:,j+1)=y(:,j)+(1/6)*(k1+2*k2+2*k3+k4);
    end
    hold on
    plot(y(1,:),y(2,:))
end
hold off