function rk4(f, t0, y0, t1, n) {
	h = (t1-t0)/(n-1)
	a = J(n, 2, 0)
	a[1, 1] = t = t0
	a[1, 2] = y = y0
	for (i=2; i<=n; i++) {
		k1 = h*(*f)(t, y)
		k2 = h*(*f)(t+0.5*h, y+0.5*k1)
		k3 = h*(*f)(t+0.5*h, y+0.5*k2)
		k4 = h*(*f)(t+h, y+k3)
		t = t+h
		y = y+(k1+2*k2+2*k3+k4)/6
		a[i, 1] = t
		a[i, 2] = y
	}
	return(a)
}

function f(t, y) {
	return(t*sqrt(y))
}

a = rk4(&f(), 0, 1, 10, 101)
t = a[., 1]
a = a, a[., 2]:-(t:^2:+4):^2:/16
a[range(1,101,10), .]

                   1              2              3
     +----------------------------------------------+
   1 |             0              1              0  |
   2 |             1    1.562499854   -1.45722e-07  |
   3 |             2    3.999999081   -9.19479e-07  |
   4 |             3    10.56249709   -2.90956e-06  |
   5 |             4    24.99999377   -6.23491e-06  |
   6 |             5    52.56248918   -.0000108197  |
   7 |             6    99.99998341   -.0000165946  |
   8 |             7    175.5624765   -.0000235177  |
   9 |             8    288.9999684   -.0000315652  |
  10 |             9    451.5624593   -.0000407232  |
  11 |            10     675.999949   -.0000509833  |
     +----------------------------------------------+