Hi everyone,

I started a thread before called "Lorenz Equations MATLAB", and this is based on the question, but im having a different problem.

So from before,

I have the following system of differential 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 < 20.

Here is my matlab code for carrying this out (contains 3 seperate mfiles)

mfile 1

----------------------------------------------------------------------

function z = ydot(t,y)

% Lorenz Equations

% z: derivative of state y

s=10; r=28; b=8/3;

z=zeros(3,1);

z(1) = s*(y(2) - y(1));

z(2) = y(1)*(r-y(3))-y(2);

z(3) = y(1)*y(2) - b*y(3);

----------------------------------------------------------------------

mfile 2

----------------------------------------------------------------------

function y = trapezoid(y,h)

z1 = ydot(0,y);

g=y+h*z1';

z2=ydot(0,g);

y=y+ 0.5*h*(z1+z2)';

---------------------------------------------------------------------

mfile 3

---------------------------------------------------------------------

function [t,y] = lorentz(int,ic,h)

%

% this function solves and plots the Lorenz equation

% t : time step vector

% y : state variable (coordinates)

%

% int: integration interval

% ic: initial state coordinates

% h: integration step size

t = 0:h:20; % set up integration time steps

n = numel(t);

y = zeros(n,3); % preallocate space

y(1, : )=ic;

t(1) = int(1);

for k=2:n

y(k, : ) = trapezoid(y(k-1, : ),h);

end

figure(1); %bring figure 1 to the front

subplot(1,2,1)

plot3(y(:,1),y(:,2 ),y(:,3));

xlabel('x');

zlabel('z');

ylabel('y');

subplot(1,2,2)

plot(y(:,1),y(:,2)),grid on;

xlabel('y1');ylabel('y2');

title('Lorenz Attractor');

-----------------------------------------------------------------------

Now to run this, in the command window type

lorentz([0 20], [5 5 5], 0.001)

What I have to do now is calculate the difference between the two trajectories and .

I thought I could do this by typing in the command window

lorentz([0 20], [5 5 5], 0.001) - lorentz([0 20], [(5 + 0.0001) 5 5], 0.001)

but this returns a value of zero.

Any ideas on how to do this? Thanks in advance.