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