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