RK4 code doesn't produce the expected output

Question

I'm trying to program my own ode4 method for solving a system of differential equations which is :

x′ = σy x(t0 = 0) = 0
y′ = ρy + βx y(t0 = 0) = 0.01

The code doesn't produce the expected output and I can't see where the error comes from. In this case, tout = tspan. I transposed my k values because the dimensions didn't match and the code didn't work without them. I have ploted the result that should be obtain with 2 other method ode45 and ode23.

Here's the code (2 files, one for the parameters and the system and one for the ode4 function).

clear all
clc



function dtdy=kernel(t,y,sigma,beta,rho)
  dtdy=zeros(2,1);
  dtdy(1)=sigma*y(2);
  dtdy(2)=rho*y(2)+beta*y(1);
end

sigma=1;
rho=-0.9;
beta=-25;
h=0.08;
t=0:h:10;
y0=[0 0.01];

[tRk4,vRk4]=ode4(@(t,y)kernel(t,y,sigma,beta,rho),t,y0);
[t23,v23]=ode23(@(t,y)kernel(t,y,sigma,beta,rho),t,y0);
[t45,v45]=ode45(@(t,y)kernel(t,y,sigma,beta,rho),t,y0);

%plot
plot(vRk4(:,1),vRk4(:,2),"b-o",v23(:,1),v23(:,2),'g-^',v45(:,1),v45(:,2),'r-')
 set(gca, "linewidth", 0.8, "fontsize", 20)
    title('Solution of a system of ODEs by Runge-Kutta methods of different order','fontsize',24);
    xlabel('x',"fontsize", 22);
    ylabel('y',"fontsize", 22);
    legend("4th order RK (ode4, my own)", "3rd order (ode23)","5th order (ode45)","fontsize",18)

function[tout,v]=ode4(odefun,tspan,v0)
v=zeros(length(tspan),length(v0));
v(1,:)=v0;%initial condition for y(0)
h=tspan(2)-tspan(1); % extract h value from tspan
tout=tspan;

for i=1:length(tspan)-1

  k1=odefun(tspan(i),v(i,:));
  k2=odefun(tspan(i)+(h/2),v(i,:)+(h/2)*k1);
  k3=odefun(tspan(i)+(h/2),v(i,:)+(h/2)*k2);
  k4=odefun(tspan(i)+h,v(i,:)+h*k3);
  v(i+1,:)=v(i,:)+(h/6)*(k1'+(2*k2')+(2*k3')+k4');
end
end

this the output : We can see than RK4 start from the good initial values but is going in the wrong way then

Answer 1

Matlab/Octave knows arrays in one dimension, and row and column vectors as objects with a 2D shape. However, these 2D vectors also can be addressed as arrays, with matrices you have to care that the Fortran column-first data arrangement is used.

Matlab/Octave and numpy in their tradition were designed to be "helpful". In certain situations a dimension mismatch is corrected via a technique called broadcasting. This allows to add row and column vectors, the input vectors just get replicated until the resulting matrices have compatible dimensions.

This is what happened here. While you corrected the situation where the incompatibility was exposed as error, you left similar mismatches where they were silently "corrected" in the broadcasting way.

To recapitulate: y(i,:) is a row vector, while your odefun produces column vectors, explicitly constructed as such with dtdy=zeros(2,1). So v(i,:)+(h/2)*k1 gets "corrected" via broadcasting as a 2x2 matrix, which gives at least one wrong value.

You could correct this at the construction or operation level, but better leave the function that is also used by the other solvers unchanged, change to

  k1=odefun(tspan(i),v(i,:));
  k2=odefun(tspan(i)+(h/2),v(i,:)+(h/2)*k1');
  k3=odefun(tspan(i)+(h/2),v(i,:)+(h/2)*k2');
  k4=odefun(tspan(i)+h,v(i,:)+h*k3');

or perhaps more systematically, that is, with less transpositions

  k1=odefun(tspan(i),v(i,:))';
  k2=odefun(tspan(i)+(h/2),v(i,:)+(h/2)*k1)';
  k3=odefun(tspan(i)+(h/2),v(i,:)+(h/2)*k2)';
  k4=odefun(tspan(i)+h,v(i,:)+h*k3)';
  v(i+1,:)=v(i,:)+(h/6)*(k1+(2*k2)+(2*k3)+k4);

This gives the expected plot of visually indistinguishable numerical solutions.