See Rosetta Code.
tStart = 0.0 tStop =10.0 dt = 0.1 ' 0.125 IS expressible exactly in binary, unlike 0.1 y = 1.0 t =tStart while t <tStop k1 = yDot( t , y ) k2 = yDot( t +0.5 *dt, y +0.5 *dt *k1) k3 = yDot( t +0.5 *dt, y +0.5 *dt *k2) k4 = yDot( t + dt , y + dt *k3) yAnalytic =( ( t^2 +4)^2) /16 if abs( int( t) -t) <dt then print "y( "; using( "##.#", t); ") ="; using( "####.#######", y), "Error ="; yAnalytic -y y = y + dt *( k1 + 2 *( k2 + k3 ) + k4 ) /6 t =t +dt wend end function yDot( t, y) yDot =t *y^0.5 end function