How to solve 4th order Runge-Kutta for different initial conditions?

I have a code that solves the 2 populations with 1 initial conditions and just plot that. How do I make 2D plots for different initial conditions?

%dx/dt = 5*x+1*x*y 
%dy/dt = -2*y+1*x*y 

clear
clc

f=@(t,z) [5*z(1)+1*z(1)*z(2);
          -2*z(2)+1*z(1)*z(2)];

% Initial Conditions
t(1)=0;
z(:,1)=[1,1];

% Step size
s=1;
e=10;
N=e/s;

% Update Loop
for i=1:N
    % Update time
    t(i+1)=t(i)+s;
    % Update for y
    k1=f(t(i)    ,z(:,i));
    k2=f(t(i)+s/2,z(:,i)+s/2*k1);
    k3=f(t(i)+s/2,z(:,i)+s/2*k2);
    k4=f(t(i)+s  ,z(:,i)+s  *k3);

    z(:,i+1)=z(:,i)+h/6*(k1 + 2*k2  + 2*k3 + k4);
end

[Davide Masiello]

Your code does not run because h is not defined.

Moreover, what exactly would you like to plot in 2D?

Answer 1

[Alan Stevens]

Here's a rather crude method (together with some corrections). You should be able to turn this into a much more elegant version:

%dx/dt = 5*x+1*x*y 
%dy/dt = -2*y+1*x*y 

f=@(t,z) [5*z(1)+1*z(1)*z(2);
          -2*z(2)+1*z(1)*z(2)];

% Example initial conditions      
x0 = [1, 0.5, 0.1];
y0 = [-1, -0.5, -0.1];


% Loop through different initial conditions
for j = 1:numel(x0)  
    % Initial Conditions
    z(:,1)=[x0(j),y0(j)];   %%%  Make y(0) negative

    % Step size
    s=0.1;           %%% Need much smaller step size
    e=10;
    N=e/s;
    
    t = zeros(1,N);
    % Update Loop
    for i=1:N
        % Update time
        t(i+1)=t(i)+s;
        % Update for y
        k1=f(t(i)    ,z(:,i));
        k2=f(t(i)+s/2,z(:,i)+s/2*k1);
        k3=f(t(i)+s/2,z(:,i)+s/2*k2);
        k4=f(t(i)+s  ,z(:,i)+s  *k3);

        z(:,i+1)=z(:,i)+s/6*(k1 + 2*k2  + 2*k3 + k4);   %%% s not h
    end

    x = z(1,:); y = z(2,:);

    figure
    subplot(2,1,1)
    plot(t,x),grid
    xlabel('t'),ylabel('x')
    subplot(2,1,2)
    plot(t,y),grid
    xlabel('t'),ylabel('y')
end

[Kenny Gee]

Thank you for answering.

I was wondering how do I plot the initial values together in one graph?

For example, I want to graph all 3 inital values for the x values under one graph. But I wanted to color each initial values differently and make a legend saying blue=initial value 1, red=initial value 2 etc. But I don't know how to separate each initial values to make a legend for them.

    figure(1); 
    
    plot(t,z(1,:),'g')
    hold on
    xlabel('Time')
    ylabel('Populations')
    legend('IV1','IV2')
    set(gca,'Fontsize',16)

[Alan Stevens]

Like this?

%dx/dt = 5*x+1*x*y 
%dy/dt = -2*y+1*x*y 

f=@(t,z) [5*z(1)+1*z(1)*z(2);
          -2*z(2)+1*z(1)*z(2)];
% Step size
s=0.1;           
e=10;
N=e/s;

% Example initial conditions      
x0 = [1, 0.5, 0.1];
y0 = [-1, -0.5, -0.1];
t = zeros(1,N);
x = zeros(numel(x0),N);   %%%%%%%%%%
y = zeros(numel(y0),N);   %%%%%%%%%%


% Loop through different initial conditions
for j = 1:numel(x0)  
        
    z(1,1) = x0(j);
    z(2,1) = y0(j);
   
    % Update Loop
    for i=1:N-1
        % Update time
        t(i+1)=t(i)+s;
        % Update for y
        k1=f(t(i)    ,z(:,i));
        k2=f(t(i)+s/2,z(:,i)+s/2*k1);
        k3=f(t(i)+s/2,z(:,i)+s/2*k2);
        k4=f(t(i)+s  ,z(:,i)+s  *k3);

        z(:,i+1)=z(:,i)+s/6*(k1 + 2*k2  + 2*k3 + k4);   %%% s not h
        
    end

    x(j,:) = z(1,:); y(j,:) = z(2,:);   %%%%%%%%%%%%%  

end
lbl1 = ['(' num2str(x0(1)) ',' num2str(y0(1)) ')'];
lbl2 = ['(' num2str(x0(2)) ',' num2str(y0(2)) ')'];
lbl3 = ['(' num2str(x0(3)) ',' num2str(y0(3)) ')'];
figure
plot(t,x),grid
xlabel('t'),ylabel('x')
legend(lbl1,lbl2,lbl3)

figure
plot(t,y),grid
xlabel('t'),ylabel('y')
legend(lbl1,lbl2,lbl3)

[Kenny Gee]

yes! thank you so much!!!