$$\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);
$(sol)-\frac{1}{y}=\frac{{{x}^{2}}}{2}+\mathit{\%c}$
 (%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.