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

 x' = 10(y - x)

 y' = x(28 - z) -y

 z' = xy - \frac{8}{3}z

I need to calculate the trajectory of  Y(t) = (x(t), y(t), z(t)) in MATLAB using Heun's method (which I think is also called the trapezoidal method) with step size h=0.001, and initial conditions  Y(0) = (5,5,5) 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  Y_1(0) = (5,5,5) and  Y_1(0) = ((5+0.0001) ,5,5) .

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.