I use MATLAB to solve the following Lorenz initial value problem:

I wrote a function, LorenzRK4IVP(), that takes the system of three differential equations as input and solves the system using the Runge-Kutta method with step size . I plot the strange attractor as well as use MATLAB to produce a GIF of the solution.

% Inputs:
% f1,f2,f3 = y'(t,y) as a string
% y0 = initial condition
% inter = interval
% h = step size
function [a] = LorenzRK4IVP(f1,f2,f3,inter,y0,h)
t(1)=inter(1);
y(1,:)=y0;
n=(inter(2)-inter(1))/h;
% evaluate derivative
function z=ydot(f1,f2,f3,t,y)
g1=inline(f1,'y1','y2','y3','t');
g2=inline(f2,'y1','y2','y3','t');
g3=inline(f3,'y1','y2','y3','t');
z(1)=g1(y(1),y(2),y(3),t(1));
z(2)=g2(y(1),y(2),y(3),t(1));
z(3)=g3(y(1),y(2),y(3),t(1));
end
% evaluate Runge-Kutta 4 formula
function y=trap(f1,f2,f3,t,y,h)
s1=ydot(f1,f2,f3,t,y);
s2=ydot(f1,f2,f3,t+(h/2),y+(h/2)*s1);
s3=ydot(f1,f2,f3,t+(h/2),y+(h/2)*s2);
s4=ydot(f1,f2,f3,t+(h),y+(h)*s3);
y=y+(h/6)*(s1+2*s2+2*s3+s4);
end
% loop
for i=1:n
t(i+1)=t(i)+h;
y(i+1,:)=trap(f1,f2,f3,t(i),y(i,:),h);
end
a=[t',y(:,1),y(:,2),y(:,3)];
end
% call function to solve Lorenz equations
L=LorenzRK4IVP('-10*y1+10*y2','-y1*y3+28*y1-y2','y1*y2-(8/3)*y3',[0,50],[5,5,5],.01)
% plot Lorenz solutions
plot3(L(:,2),L(:,3),L(:,4))
% create GIF source:
% http://www.mathworks.com/matlabcentral/answers/94495-how-can-i-create-animated-gif-images-in-matlab
figure(1)
filename = 'Lorenz.gif';
for n = 1:(length(L)/10)
plot3(L(1:10*n,2),L(1:10*n,3),L(1:10*n,4))
drawnow
frame = getframe(1);
im = frame2im(frame);
[imind,cm] = rgb2ind(im,256);
if n == 1;
imwrite(imind,cm,filename,'gif', 'Loopcount',inf);
else
imwrite(imind,cm,filename,'gif','WriteMode','append');
end
end
L=LorenzRK4IVP('-10*y1+10*y2','-y1*y3+28*y1-y2','y1*y2-(8/3)*y3',[0,50],[5,5,5],.01)