G'day all,

I have the following Lorenz equations

I need to calculate the trajectory of in MATLAB using Heun's method (which I think is also called the trapezoidal method) with step size h=0.001, and initial conditions for t < 30.

Here is my matlab code for a single IVP (using Heun's method):

function y = heuns(f,t,x,h)

g = f(t,x);

z = x + h * g;

y = x + 0.5 * h * ( g + f(t+h,z) );

where f is a MATLAB inline function, h is the step size, t is the initial time, and x is the initial condition. However, I dont know how to use this method for the system of Lorenz equations above. Any suggestions?