first test the convergence in space by taking dt fixed and making dx
smaller by halves:
>> [u,err1,x1,t1] = heat1(0,.0001,100,100);
r =
  2.4826e-004

>> [u,err2,x2,t2] = heat1(0,.0001,100,200);
r =
    0.0010

I have three solutions and three errors at a bunch of mesh-points.  I
want to compare errors at a mesh-point that is in all three sets of
mesh-points.  The first mesh-point, x=0, is in all three sets.  So I
compare the errors at x=0 over the entire time window.  The space
derivatives are O(dx^2) so the ratios should be close to 4 if the
code is right.
>> rat = err1(1,:)./err2(1,:);
Warning: Divide by zero.

>> format long 
>> min(rat), max(rat)
ans =
   4.05850507060528
ans =
   4.05850655281055

I can look at the maximum error over the entire space-time grid
>> max(max(abs(err1)))
ans =
    3.350883370778490e-008

>> max(max(abs(err2)))
ans =
    8.256444972154498e-009

the ratio of the maximum errors is also close to 4.  This is an
illegitimate test, of course since the maximum errors could be
occurring at different space-time mesh-points and hence shouldn't
be compared to each other.
>> max(max(abs(err1)))/max(max(abs(err2)))
ans =
   4.05850627246909

The above demonstrates that the space derivatives have been done
correctly.  

Now we check the time-derivatives.  This is a little tricky.  First, I
take a fairly fine space mesh.  I then take the time-steps smaller and
smaller until I find that the error bottoms out around 1.e-8.  That
is, I can't make the error smaller by making the time-steps smaller.
To make it smaller, I'd have to take the space-mesh finer.  Now that I
know that it bottoms out around 1.e-8, I choose the time steps large
enough that the errors are significantly larger than this so the
errors can decrease as they should.

>> [u,err1,x,t] = heat1(0,.1,1,64);
r =
  1.0054e+001
>> max(max(abs(err1)))
ans =
  4.7546e-003

>> [u,err2,x,t] = heat1(0,.1,2,64);
r =
  5.0268e+000
>> max(max(abs(err2)))
ans =
  2.2587e-003

>> [u,err3,x,t] = heat1(0,.1,4,64);
r =
  2.5134e+000
>> max(max(abs(err3)))
ans =
  1.0727e-003

>> [u,err4,x,t] = heat1(0,.1,8,64);
r =
  1.2567e+000

Now I test the ratios at the final time, t=.1 They should be close to
2, since the time-stepping is O(dt).
>> max(err1(:,1+1)./err2(:,2+1)) , min(err1(:,1+1)./err2(:,2+1)) 
ans =
  2.1050e+000
ans =
  2.1050e+000
>> max(err2(:,2+1)./err3(:,4+1)) , min(err2(:,2+1)./err3(:,4+1)) 
ans =
  2.1056e+000
ans =
  2.1056e+000
>> max(err3(:,4+1)./err4(:,8+1)) , min(err3(:,4+1)./err4(:,8+1)) 
ans =
  2.1705e+000
ans =
  2.1705e+000

This confirms that the time-stepping is O(dt) and the spatial derivatives
are done with O(dx^2) accuracy.   
>> diary off


%
>> [u1,err1,x,t] = heat1(0,.001,10,500);
r =
    0.6307
>> [u2,err2,x,t] = heat1(0,.001,20,500);
r =
    0.3154
>> [u3,err3,x,t] = heat1(0,.001,40,500);
r =
    0.1577
>> [u4,err4,x,t] = heat1(0,.001,80,500);
r =
    0.0788
>> [u5,err5,x,t] = heat1(0,.001,160,500);
r =
    0.0394
>> [u6,err6,x,t] = heat1(0,.001,320,500);
r =
    0.0197
>> rat1 = (u1(:,11)-u2(:,21))./(u2(:,21)-u3(:,41));
>> rat2 = (u2(:,21)-u3(:,41))./(u3(:,41)-u4(:,81));
>> max(rat1), min(rat1)
ans =
    2.0002
ans =
    2.0001
>> max(rat2), min(rat2)
ans =
    2.0001
ans =
    2.0000

% this confirms that the time-stepping is O(dt) using the computed
% solutions.  Now we verify this using the errors.
>> [u6,err6,x,t] = heat1(0,.0001,320,320);
r =
   8.0551e-04
>> [u7,err7,x,t] = heat1(0,.0001,640,320);
r =
   4.0276e-04
>> rat = err6(:,321)./err7(:,321);
>> plot(rat)
>> max(rat), min(rat)
ans =
    1.9951
ans =
    1.9950


>> [u1,err1,x,t] = heat1(0,.0001,320,10);
r =
   6.4117e-07
>> [u2,err2,x,t] = heat1(0,.0001,320,20);
r =
   2.8576e-06
>> [u3,err3,x,t] = heat1(0,.0001,320,40);
r =
   1.2040e-05
>> [u4,err4,x,t] = heat1(0,.0001,320,80);
r =
   4.9402e-05
>> [u5,err5,x,t] = heat1(0,.0001,320,160);
r =
   2.0012e-04
>> [u6,err6,x,t] = heat1(0,.0001,320,320);
r =
   8.0551e-04
>> rat1 = (u1(1,:)-u2(1,:))./(u2(1,:)-u3(1,:));
Warning: Divide by zero.
>> rat2 = (u2(1,:)-u3(1,:))./(u3(1,:)-u4(1,:));
Warning: Divide by zero.

>> rat3 = (u3(1,:)-u4(1,:))./(u4(1,:)-u5(1,:));
Warning: Divide by zero.
>> rat4 = (u4(1,:)-u5(1,:))./(u5(1,:)-u6(1,:));
Warning: Divide by zero.
>> max(rat1), min(rat1)
ans =
    4.4633
ans =
    4.4633
>> max(rat2), min(rat2)
ans =
    4.2342
ans =
    4.2342
>> max(rat3), min(rat3)
ans =
    4.1171
ans =
    4.1171
>> max(rat4), min(rat4)
ans =
    4.0585
ans =
    4.0584

% this confirms the spatial discretization is O(dx^2), using the solution
% values.  Now we test with the errors:

>> rat = err5(1,:)./err6(1,:);
Warning: Divide by zero.
>> max(rat), min(rat)
ans =
    4.0397
ans =
    4.0397

>> diary off