Runga-Kutta Method for system of first order differential equations

I am trying to work out the runga kutta method for a system of first order ODE's. I have the scheme already set up and it runs without any errors, giving an estimation. The ODE's I am trying to solve:

This scheme blows up, which is not what I was expecting. Is there an error in the scheme? Or have I missunderstood the fact that this scheme should not blow up?

Below you can find the code I have written:

%% Initialization

tspan = 100;

h = 0.125;

fact = 1/h;

x = 0:h:tspan;

y = zeros(length(x),2);

y(1,1) = 0.2; % Initial condition y1(0) = 0.2

y(1,2) = 0; % Initial coniditon y2(0) = 0

dy1dt = @(t,y) y(:,2); % change the function as you desire

dy2dt = @(t,y) -y(:,1); % change the function as you desire

%% Runge Kutta

if Method == 'RKM' | Method == 'All'

for i3 = 1:(length(x)-1) % calculation loop

k_1 = dy1dt(x(i3) , y(i3,:));

k_2 = dy1dt(x(i3) + 0.5*h , y(i3,:) - 0.5.*h.*k_1);

k_3 = dy1dt((x(i3) + 0.5*h) , (y(i3,:) - 0.5.*h.*k_2));

k_4 = dy1dt((x(i3) + h) , (y(i3,:) - k_3.*h));

y(i3+1,1) = y(i3,1) + (1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4).*h; % main equation

l_1 = dy2dt(x(i3) , y(i3,:));

l_2 = dy2dt(x(i3) + 0.5*h , y(i3,:) - 0.5*h*l_1);

l_3 = dy2dt((x(i3) + 0.5*h) , (y(i3,:) - 0.5*h*l_2));

l_4 = dy2dt((x(i3) + h) , (y(i3,:) - l_3*h));

y(i3+1,2) = y(i3,2) + (1/6)*(l_1 + 2*l_2* + 2*l_3 + l_4)*h; % main equation

end

figure(3),clf(3),hold on

plot(x,y(:,1),'r')

plot(x,y(:,2),'b')

plot(T,Y(:,1),'r--') % Solution according to ODE45

plot(T,Y(:,2),'b--') % Solution according to ODE45

axis([0 100 -1 1])

legend('Runga Kutta 1','Runga Kutta 2','Exact 1','Exact 2')

title('Runga Kutta')

grid on;

end

The result of this code is below.

Answer 1

[Davide Masiello]

Two things:

1 - By defining two different functions for y1 and y2 you decouple the system, which is the source of the main problem.

2 - After fixing point 1, I had to use a very small time step size to obtain an acceptable solution, indicating that the implementation of an adaptive step size might be beneficial to the efficiency of the code.

%% Initialization
tspan   = 100;
h       = 0.0005;
fact    = 1/h;
x       = 0:h:tspan;
y       = zeros(2,length(x));
y(:,1)  = [0.2;0];                     % Initial conditions

%% Runge Kutta
for i = 1:(length(x)-1)                              % calculation loop
    k_1         = odeSystem(x(i),y(:,i));
    k_2         = odeSystem(x(i)+0.5*h,y(:,i)-0.5.*h.*k_1);
    k_3         = odeSystem((x(i)+0.5*h),(y(:,i)-0.5*h*k_2));
    k_4         = odeSystem((x(i)+h),(y(:,i)-k_3*h));
    y(:,i+1)    = y(:,i)+(1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4)*h;  % main equation
end
    
figure
plot(x,y(1,:),'r',x,y(2,:),'b')
axis([0 100 -1 1])
title('Runge-Kutta')

function dydt = odeSystem(t,y)
dydt(1,1) = y(2); 
dydt(2,1) = -y(1);
end

Answer 2

[Chunru]

Looks like you need to choose a much smaller step size for the problem.

%% Initialization
tspan = 100;
h = 0.125/50;
fact = 1/h;
x = 0:h:tspan;
y = zeros(length(x),2);
y(1,1) = 0.2;                     % Initial condition y1(0) = 0.2
y(1,2) = 0;                       % Initial coniditon y2(0) = 0
dy1dt = @(t,y) y(:,2);                    % change the function as you desire
dy2dt = @(t,y) -y(:,1);                    % change the function as you desire
%% Runge Kutta
Method = "RKM"

Method = "RKM"

if Method == 'RKM' | Method == 'All'
    for i3 = 1:(length(x)-1)                              % calculation loop
        k_1 = dy1dt(x(i3) , y(i3,:));
        k_2 = dy1dt(x(i3) + 0.5*h , y(i3,:) - 0.5.*h.*k_1);
        k_3 = dy1dt((x(i3) + 0.5*h) , (y(i3,:) - 0.5.*h.*k_2));
        k_4 = dy1dt((x(i3) + h) , (y(i3,:) - k_3.*h));
        y(i3+1,1) = y(i3,1) + (1/6)*(k_1 + 2*k_2 + 2*k_3 + k_4).*h;  % main equation
        l_1 = dy2dt(x(i3) , y(i3,:));
        l_2 = dy2dt(x(i3) + 0.5*h , y(i3,:) - 0.5*h*l_1);
        l_3 = dy2dt((x(i3) + 0.5*h) , (y(i3,:) - 0.5*h*l_2));
        l_4 = dy2dt((x(i3) + h) , (y(i3,:) - l_3*h));
        y(i3+1,2) = y(i3,2) + (1/6)*(l_1 + 2*l_2* + 2*l_3 + l_4)*h;  % main equation
    end
    
    figure(3),clf(3),hold on
    plot(x,y(:,1),'r')
    plot(x,y(:,2),'b')
    
%     plot(T,Y(:,1),'r--')       % Solution according to ODE45
%     plot(T,Y(:,2),'b--')       % Solution according to ODE45
    
    axis([0 100 -1 1])
    legend('Runga Kutta 1','Runga Kutta 2','Exact 1','Exact 2')
    title('Runga Kutta')
    grid on;
end

Warning: Ignoring extra legend entries.