function lorenz_rk

%-----------------

N=5000;

sigma = 10;

rho = 25;

beta = 8/3.;

t=linspace(0,50,N); % time between 0-5s

h=t(2)-t(1);   % time step

f=@(x,y) sigma*(y-x);

g=@(x,y,z) x*(rho-z)-y;

u=@(x,y,z) x*y-beta*z;

x(1)=1;

y(1)=0;

z(1)=0;

for i=1:length(t)-1

kx1=f(x(i),y(i));

ky1=g(x(i),y(i),z(i));

kz1=u(x(i),y(i),z(i));

%------------------------

kx2=f(x(i)+h/2*kx1,y(i)+h/2*ky1);

ky2=g(x(i)+h/2*kx1,y(i)+h/2*ky1,z(i)+h/2*kz1);

kz2=u(x(i)+h/2*kx1,y(i)+h/2*ky1,z(i)+h/2*kz1);

%------------------------

kx3=f(x(i)+h/2*kx2,y(i)+h/2*ky2);

ky3=g(x(i)+h/2*kx2,y(i)+h/2*ky2,z(i)+h/2*kz2);

kz3=u(x(i)+h/2*kx2,y(i)+h/2*ky2,z(i)+h/2*kz2);

%------------------------

kx4=f(x(i)+h*kx3,y(i)+h*ky3);

ky4=g(x(i)+h*kx3,y(i)+h*ky3,z(i)+h*kz3);

kz4=u(x(i)+h*kx3,y(i)+h*ky3,z(i)+h*kz3);

%------------------------

x(i+1)=x(i)+h/6*(kx1+2*kx2+2*kx3+kx4);

y(i+1)=y(i)+h/6*(ky1+2*ky2+2*ky3+ky4);

z(i+1)=z(i)+h/6*(kz1+2*kz2+2*kz3+kz4);

end

body = line( ...

'color','b', ...

'LineStyle','none', ...

'Marker','.', ...

'erase','none', ...

'xdata',[],'ydata',[],'zdata',[]);

axis([ min(y)-eps max(y)+eps...

min(z)-eps max(z)+eps])

for j=2:N

set(body,'XData',y(j),'YData',z(j))

drawnow

end

plot(y,z,'.')

xlabel('y')

ylabel('z')

function Lorenz_ode

close all

x0=1;

y0=0;

z0=0;

N=10000;

[t,Y] = ode45(@lorenz_at,linspace(0,50,N),[x0 y0 z0]);

body = line( ...

'color','b', ...

'LineStyle','none', ...

'Marker','.', ...

'erase','none', ...

'xdata',[],'ydata',[],'zdata',[]);

axis([ min(Y(:,2))-eps max(Y(:,2))+eps...

min(Y(:,3))-eps max(Y(:,3))+eps])

for j=2:N

set(body,'XData',Y(j,2),'YData',Y(j,3))

drawnow

end

plot(Y(:,2),Y(:,3),'.')

xlabel('y')

ylabel('z')

function dy=lorenz_at(t,y)

dy=zeros(3,1);

SIGMA = 10;

RHO = 25;

BETA = 8/3.;

dy(1)=SIGMA*(y(2)-y(1));

dy(2)=y(1)*(RHO-y(3))-y(2);

dy(3)=y(1)*y(2)-BETA*y(3);