MATLAB: Solving Lorenz attractor equations using Runge kutta (RK4) method

MATLAB

I am trying to write a code for the simulation of lorenz attractor using rk4 method. Here is the code:
clc;
clear all;
t(1)=0;  %initializing x,y,z,t
x(1)=1;
y(1)=1;
z(1)=1;
sigma=10;   %value of constants
rho=28;
beta=(8/3);
h=0.1;   %step size
t=0:h:100;
f=@(t,x,y,z) sigma*(y-x);  %ode 
g=@(t,x,y,z) x*rho-x.*z-y;
p=@(t,x,y,z) x.*y-beta*z;
for i=1:(length(t)-1) %loop
    k1=h*f(t(i),x(i),y(i),z(i));
    l1=h*g(t(i),x(i),y(i),z(i));
    m1=h*p(t(i),x(i),y(i),z(i));
      k2=h*f(t(i)+h,(x(i)+0.5*k1),(y(i)+(0.5*l1)),(z(i)+(0.5*m1)));
      l2=h*f(t(i)+h,(x(i)+0.5*k1),(y(i)+(0.5*l1)),(z(i)+(0.5*m1)));
      m2=h*f(t(i)+h,(x(i)+0.5*k1),(y(i)+(0.5*l1)),(z(i)+(0.5*m1)));
      k3=h*f(t(i)+h,(x(i)+0.5*k2),(y(i)+(0.5*l2)),(z(i)+(0.5*m2)));
      l3=h*f(t(i)+h,(x(i)+0.5*k2),(y(i)+(0.5*l2)),(z(i)+(0.5*m2)));
      m3=h*f(t(i)+h,(x(i)+0.5*k2),(y(i)+(0.5*l2)),(z(i)+(0.5*m2)));
      k4=h*f(t(i)+h,(x(i)+k3),(y(i)+l3),(z(i)+m3));
      l4=h*g(t(i)+h,(x(i)+k3),(y(i)+l3),(z(i)+m3));
      m4=h*p(t(i)+h,(x(i)+k3),(y(i)+l3),(z(i)+m3));
      x(i+1)=x(i)+h*(1/6)*(k1+2*k2+2*k3+k4); %final equations
      y(i+1)=y(i)+h*(1/6)*(k1+2*k2+2*k3+k4);
      z(i+1)=z(i)+h*(1/6)*(m1+2*m2+2*m3+m4);
    end
 plot3(x,y,z)
But the solutions are not right. I don't know what to do. I know we can do using ode solvers but i wanted to do using rk4 method. I searched for the solutions in different sites but i didn't find many using rk4. While there were some but only algorithm. I tried to compare my solutions with ode45 but doesn't match at all. it's totally different.

Best Answer

The only bug that I can see at first glance is here
 y(i+1) = y(i) + h*(1/6)*(l1+2*l2+2*l3+l4); % replace ki by li
You also might want to play with (decrease) the step size.
Related Question