\( \DeclareMathOperator{\abs}{abs} \)

Let's check that the built-in Runge Kutta numerical ode solver is really 4th order accurate.  We'll use the nonlinear 1st order IVP \[ y'=xy^2, y(0)=\frac{1}{2}.\]
First the analytical solution for \( y(1) \):

(%i1) sol: ode2('diff(y,x)=x*y^2,y,x);
(%i2) ic1(%, x=0, y=1/2);
\[\mathrm{\tt (\%o2) }\quad -\frac{1}{y}=\frac{{{x}^{2}}-4}{2}\]
(%i3) subst(x=1,%);
\[\mathrm{\tt (\%o3) }\quad -\frac{1}{y}=-\frac{3}{2}\]
(%i4) solve(%,y);
\[\mathrm{\tt (\%o4) }\quad [y=\frac{2}{3}]\]

Now we'll run rk on the interval [0,1] with timesteps that decrease in size by factors of 10

(%i5) s1: rk(x*y^2,y,.5,[x,0,1,.1])$
(%i6) s01: rk(x*y^2,y,.5,[x,0,1,.01])$
(%i7) s001: rk(x*y^2,y,.5,[x,0,1,.001])$

An finally, compute the error in our approximation of \( y(1) \). Notice that when the stepsize is decreased by a factor of 10, the error decreases by a factor of \( 10^4 \)...fourth order accurate!

(%i8) 2/3-s1[11][2];
\[\mathrm{\tt (\%o8) }\quad -1.54354535597534\cdot {{10}^{-7}}\]
(%i9) 2/3-s01[101][2];
\[\mathrm{\tt (\%o9) }\quad -2.052613634617728\cdot {{10}^{-11}}\]
(%i10) 2/3-s001[1001][2];
\[\mathrm{\tt (\%o10) }\quad -3.441691376337985\cdot {{10}^{-15}}\]
Created with wxMaxima. The source of this maxima session can be downloaded here.