compute with time step h=1
>> [x,y] = euler_step(1,0,10,1);
check that we have the right initial time:
>> x(1)
>> ans =
     0
check that we have the right initial value:
>> y(1)
ans =
     1

define the true solution:
>> Y = 1*exp(-2*x);
plot the true solution and the computed solution on the same graph
>> plot(x,y,'x'); hold on; plot(x,Y); hold off; figure(1) 

compute with time step h=1/2
>> [x,y] = euler_step(1,0,10,1/2);
>> Y = 1*exp(-2*x);
>> size(y)
ans =
     1    21
define the error at the final time (x=10) using h=1/2:
>> err(1) = Y(21)-y(21);
plot the true solution and the computed solution on the same graph
>> plot(x,y,'x'); hold on; plot(x,Y); hold off; figure(1)

compute with time step h=1/4
>> [x,y] = euler_step(1,0,10,1/4);
>> Y = 1*exp(-2*x);
>> size(y)
ans =
     1    41
define the error at the final time (x=10) using h=1/4:
>> err(2) = Y(41)-y(41);
plot the true solution and the computed solution on the same graph
>> plot(x,y,'x'); hold on; plot(x,Y); hold off; figure(1)

compute with time step h=1/8
>> [x,y] = euler_step(1,0,10,1/8);
>> Y = 1*exp(-2*x);
>> size(y)
ans =
     1    81
define the error at the final time (x=10) using h=1/8:
>> err(3) = Y(81)-y(81);
plot the true solution and the computed solution on the same graph
>> plot(x,y,'x'); hold on; plot(x,Y); hold off; figure(1)

The computed solutions are clearly converging to the true solution
on the graphs, so I'll stop plotting them and will only define some
more errors to work with.

>> [x,y] = euler_step(1,0,10,1/16);
>> size(y)
ans =
     1    161
>> Y = 1*exp(-2*x);
>> err(4) = Y(161)-y(161);

>> [x,y] = euler_step(1,0,10,1/32);
>> size(y)
ans =
     1    321
>> Y = 1*exp(-2*x);
>> err(5) = Y(321)-y(321);

>> [x,y] = euler_step(1,0,10,1/64);
>> size(y)
ans =
     1    641
>> Y = 1*exp(-2*x);
>> err(6) = Y(641)-y(641);

>> [x,y] = euler_step(1,0,10,1/128);
>> size(y)
ans =
     1    1281
>> Y = 1*exp(-2*x);
>> err(7) = Y(1281)-y(1281);


now look at the errors:
>> err
err =
  Columns 1 through 6 
  2.0612e-009  2.0602e-009  1.9600e-009  1.5348e-009  9.8764e-010  5.6320e-010
  Column 7 
  3.0106e-010

They're decreasing to zero.  Use them to define the ratios:
>> for i=1:6, rat(i) = err(i)/err(i+1); end
>> rat
rat =
  1.0004e+000  1.0511e+000  1.2771e+000  1.5540e+000  1.7536e+000  1.8707e+000

The ratio appears to be increasing, but to what?  Let's try computing
with smaller time-steps.  But I'm too impatient to compute all the way
out to the final time of 1.

>> [x,y] = euler_step(1,0,1,1/128);
>> size(y)
ans =
     1    129
>> Y = 1*exp(-2*x);
>> clear err
>> err(1) = Y(129)-y(129);
>> y1 = y(129);

>> [x,y] = euler_step(1,0,1,1/256);
>> size(y)
ans =
     1    257
>> Y = 1*exp(-2*x);
>> err(2) = Y(257)-y(257);
>> y2 = y(257);

>> [x,y] = euler_step(1,0,1,1/512);
>> size(y)
ans =
     1    513
>> Y = 1*exp(-2*x);
>> err(3) = Y(513)-y(513);
>> y3 = y(513);

>> err
err =
  2.1201e-003  1.0587e-003  5.2900e-004
The error is decreasing.

>> err(1)/err(2)
ans =
  2.0026e+000
The first ratio is close to 2

>> err(2)/err(3)
ans =
  2.0013e+000
The second ratio is even closer to 2.

I saved the three approximate values at the final time x=1, so I could
use them to make the ratio that doesn't require knowing the true
solution:

>> y1, y2, y3
y1 =
  1.3322e-001
y2 =
  1.3428e-001
y3 =
  1.3481e-001

The values appear to be converging.

>> (y1-y2)/(y2-y3)
ans =
  2.0039e+000

The ratio is close to 2.  So the evidence shows that Euler's method
has an error that decreases like C h, rather than C h^2 or C h^3...

>> diary off