How to test your convergence graphically? The following is for dy/dx = cos(x)*y I'll compute the solutions from x=1 to x=20 with initial data y(1) = 2. I'll compute six solutions where each has h divided by 2. >> y_0 = 2; x_0 = 1; x_final = 20; N = 100; >> [x1,y1] = euler_step(y_0,x_0,x_final,N,'f'); >> N = 2*N; >> [x2,y2] = euler_step(y_0,x_0,x_final,N,'f'); >> N = 2*N; >> [x3,y3] = euler_step(y_0,x_0,x_final,N,'f'); >> N = 2*N; >> [x4,y4] = euler_step(y_0,x_0,x_final,N,'f'); >> N = 2*N; >> [x5,y5] = euler_step(y_0,x_0,x_final,N,'f'); >> N = 2*N; >> [x6,y6] = euler_step(y_0,x_0,x_final,N,'f'); Looking at the solutions, I see that the dimensions are reasonable... >> whos x* y* Name Size Bytes Class x1 1x101 808 double array x2 1x201 1608 double array x3 1x401 3208 double array x4 1x801 6408 double array x5 1x1601 12808 double array x6 1x3201 25608 double array x_0 1x1 8 double array x_final 1x1 8 double array y1 1x101 808 double array y2 1x201 1608 double array y3 1x401 3208 double array y4 1x801 6408 double array y5 1x1601 12808 double array y6 1x3201 25608 double array y_0 1x1 8 double array Grand total is 12838 elements using 102668 bytes Check that I didn't make any mistakes by checking that I did start at x=1 and stop at x=20: >> max(x1), min(x1) ans = 20 ans = 1 >> max(x2), min(x2) ans = 20 ans = 1 >> max(x3), min(x3) ans = 20 ans = 1 I know that x1 goes by steps of size h, x2 goes by steps half as long, and x3 goes by steps a quarter as long. So I can look at each step of x1, every other step of x2, and every every other step of x3. Check that the dimensions make sense: >> size(x1) ans = 1 101 >> size(x2(1:2:length(x2))) ans = 1 101 >> size(x3(1:4:length(x3))) ans = 1 101 An eyeball check to see that the errors are decreasing in magnitude: >> figure(1) >> clf >> plot(x1,abs(y1-y2(1:2:length(y2))),'b') >> hold on >> plot(x2,abs(y2-y3(1:2:length(y3))),'k') >> plot(x3,abs(y3-y4(1:2:length(y4))),'c') >> plot(x4,abs(y4-y5(1:2:length(y5))),'g') >> plot(x5,abs(y5-y6(1:2:length(y6))),'r') That's promising, but inexact. Since I know how Euler stepping is supposed to converge, I want to look at the convergence. I'll do this using as much data as I can. I would like to look at y(x_0 + h) which is approximately y1(2), y2(3), and y4(5). Since I don't know the exact solution, I will compute my convergence ratios by looking at (y1(2)-y2(3))/(y2(3)-y3(5)). This is just for x_0 + h. I can also look at this ratio at x_0 + 2*h. This would be the ratio (y1(3)-y2(5))/(y2(5)-y3(9)) and so on. >> rat1 = (y1 - y2(1:2:length(y2)))./(y2(1:2:length(y2))-y3(1:4:length(y3))); Warning: Divide by zero. (Type "warning off MATLAB:divideByZero" to suppress this warning.) >> rat2 = (y2 - y3(1:2:length(y3)))./(y3(1:2:length(y3))-y4(1:4:length(y4))); Warning: Divide by zero. (Type "warning off MATLAB:divideByZero" to suppress this warning.) >> rat3 = (y3 - y4(1:2:length(y4)))./(y4(1:2:length(y4))-y5(1:4:length(y5))); Warning: Divide by zero. (Type "warning off MATLAB:divideByZero" to suppress this warning.) >> rat4 = (y4 - y5(1:2:length(y5)))./(y5(1:2:length(y5))-y6(1:4:length(y6))); Warning: Divide by zero. (Type "warning off MATLAB:divideByZero" to suppress this warning.) What's the cause of the "divide by zero"? Well they all agree at the first point, so rat1(1) = (y_0-y_0)/(y_0-y_0). >> whos rat* Name Size Bytes Class rat1 1x101 808 double array rat2 1x201 1608 double array rat3 1x401 3208 double array rat4 1x801 6408 double array Grand total is 1504 elements using 12032 bytes Now I plot those ratios as a function of time. >> figure(2) >> clf >> plot(x1,rat1,'b') >> hold on >> plot(x2,rat2,'k') >> plot(x3,rat3,'c') >> plot(x4,rat4,'g') >> axis([x_0,x_final,0,4]) These plots are, by and large, getting closer and closer to the value 2. Which is good. Now, in fact, I actually know the exact solution for this problem. It's 2/exp(sin(1)) exp(sin(x)). When I know the exact solution, I can do something stronger than the above to test the convergence: First, I define the exact solutions >> Y1 = 2*exp(sin(x1))/exp(sin(1)); >> Y2 = 2*exp(sin(x2))/exp(sin(1)); >> Y3 = 2*exp(sin(x3))/exp(sin(1)); >> Y4 = 2*exp(sin(x4))/exp(sin(1)); >> Y5 = 2*exp(sin(x5))/exp(sin(1)); >> Y6 = 2*exp(sin(x6))/exp(sin(1)); Second I define the ratios of the errors: >> Rat1 = (y1 - Y1)./(y2(1:2:length(y2))-Y2(1:2:length(Y2))); Warning: Divide by zero. (Type "warning off MATLAB:divideByZero" to suppress this warning.) >> Rat2 = (y2 - Y2)./(y3(1:2:length(y3))-Y3(1:2:length(Y3))); Warning: Divide by zero. (Type "warning off MATLAB:divideByZero" to suppress this warning.) >> Rat3 = (y3 - Y3)./(y4(1:2:length(y4))-Y4(1:2:length(Y4))); Warning: Divide by zero. (Type "warning off MATLAB:divideByZero" to suppress this warning.) >> Rat4 = (y4 - Y4)./(y5(1:2:length(y5))-Y5(1:2:length(Y5))); Warning: Divide by zero. (Type "warning off MATLAB:divideByZero" to suppress this warning.) >> Rat5 = (y5 - Y5)./(y6(1:2:length(y6))-Y6(1:2:length(Y6))); Warning: Divide by zero. (Type "warning off MATLAB:divideByZero" to suppress this warning.) >> whos R* Name Size Bytes Class Rat1 1x101 808 double array Rat2 1x201 1608 double array Rat3 1x401 3208 double array Rat4 1x801 6408 double array Rat5 1x1601 12808 double array Grand total is 3105 elements using 24840 bytes Now I plot those ratios as a function of time. >> figure(3) >> clf >> plot(x1,Rat1,'b') >> hold on >> plot(x2,Rat2,'k') >> plot(x3,Rat3,'c') >> plot(x4,Rat4,'g') >> plot(x5,Rat5,'r') >> axis([x_0,x_final,0,4]) Why is this a stronger test? Not only does it check that it's converging at the correct *rate* it checks that the approximate solutions are converging to the correct *solution*. This is an important test to do whenever you have an exact solution handy.