1/*  Runge-Kutta Method  */
2
3rungekutta(p1, p2, p, q, tt) :=
4  block( [k11, k12, k21, k22, k31, k32, k41, k42],
5         k11 : hh*p1,
6	  k12 : hh*p2,
7	  k21 : hh*ratsubst(q+k12/2, q,
8			    ratsubst(p+k11/2, p,
9				     ratsubst(tt+hh/2, tt, p1))),
10	  k22 : hh*ratsubst(q+k12/2, q,
11			    ratsubst(p+k11/2, p,
12			 	     ratsubst(tt+hh/2, tt, p2))),
13	  k31 : hh*ratsubst(q+k22/2, q,
14			    ratsubst(p+k21/2, p,
15				     ratsubst(tt+hh/2, tt, p1))),
16	  k32 : hh*ratsubst(q+k22/2, q,
17			    ratsubst(p+k21/2, p,
18				     ratsubst(tt+hh/2, tt, p2))),
19	  k41 : hh*ratsubst(q+k32, q,
20			    ratsubst(p+k31, p,
21				     ratsubst(tt+hh, tt, p1))),
22	  k42 : hh*ratsubst(q+k32, q,
23			    ratsubst(p+k31, p,
24				     ratsubst(tt+hh, tt, p2))),
25	  pn : ratsimp(p + (k11 + 2*k21 + 2*k31 + k41)/6),
26	  qn : ratsimp(q + (k12 + 2*k22 + 2*k32 + k42)/6)
27       )$
28