--RungeKutta routine which solves simultaneous differential
--equations using the forth order Runge-Kutta approximation.
--Results appear at very bottom of this Pad

--by Anthony Egan
--SD School of Mines
--afe8607@silver.sdsmt.edu

--improved by U.Wienands <uli@SLAC.Stanford.EDU> 19-Jan-97

iters=100 --the number of divisions between x_initial and x_final
ny=3      --the number of equations

--input an array of ny differential equations.

ydot(x,y)={y[1],
          y[1]-y[2],
         x^2-x}

zeros2[i,j] = 0 dim[ny,iters+1] -- Note: iters+1 points
xx:=0[:iters+1]:
y:= 0[:ny]:
C1:=y:;C2:=y:;C3:=y:;C4:=y:
y_res:=zeros2:

--input initial conditions of the dependent variables

y[1]:=1:
y[2]:=4:
y[3]:=1:

--input the initial and final values for independent variable x

xinitial:=2:
xfinal:=6:

--Begin Program

dxx:=(xfinal-xinitial)/iters:; dxx:40.0e-3
xx[1]:=xinitial:

i:=1:
k:=1,(y_res[k,1]:=y[k],k:=k+1) while k<=ny:

(k:=1, (C1[k]:=ydot(xx[i],y)[k],k:=k+1) while k<=ny,
 k:=1, (C2[k]:=ydot(xx[i]+0.5*dxx,y+0.5*C1*dxx)[k],k:=k+1) while k<=ny,
 k:=1, (C3[k]:=ydot(xx[i]+0.5*dxx,y+0.5*C2*dxx)[k],k:=k+1) while k<=ny,
 k:=1, (C4[k]:=ydot(xx[i]+dxx,y+C3*dxx)[k],k:=k+1) while k<=ny,
 k:=1, (y_res[k,i+1]:=y[k]+dxx*(C1[k]+2*C2[k]+2*C3[k]+C4[k])/6,k:=k+1) while k<=ny,
 k:=1, (y[k]:=y_res[k,i+1],k:=k+1) while k<=ny,
xx[i+1]:=xx[i]+dxx,i:=i+1) while i<=iters:

--Begin Results

--final values

xinitial+dxx*(iters):6.00--final x
y_res[1,iters+1]:54.6--final y[1]
y_res[2,iters+1]:27.4--final y[2]
y_res[3,iters+1]:54.3--final y[3]

--plot of results.

Title:="y's vs x":
Xmin:=xinitial:;  Xmax:=xfinal:
         Ylabel:=" y[1]":;  plotline {xx,y_res[1]}
newaxis;
         Ylabel:=" y[2]":;  plotline {xx,y_res[2]}
newaxis; Ylabel:=" y[3]":;  plotline {xx,y_res[3]}

--table of other values
--if you wish to see everything, change increments to 1
increments:=5:
Nmin:=1:
Nmax:=iters+1:
Nsteps:=iters/increments:
table xx[N],y_res[1,N], y_res[2,N], y_res[3,N]
       2.00      1.000       4.00      1.000
       2.20       1.22       3.48       1.46
       2.44       1.55       3.03       2.20
       2.64       1.90       2.79       2.98
       2.84       2.32       2.67       3.94
       3.04       2.83       2.65       5.08
       3.28       3.60       2.77       6.72
       3.48       4.39       2.99       8.33
       3.68       5.37       3.34       10.2
       3.88       6.55       3.81       12.3
       4.12       8.33       4.59       15.2
       4.32       10.2       5.43       17.9
       4.52       12.4       6.50       20.9
       4.72       15.2       7.82       24.2
       4.96       19.3       9.83       28.7
       5.16       23.6       11.9       32.8
       5.36       28.8       14.5       37.3
       5.56       35.2       17.7       42.2
       5.80       44.7       22.4       48.6
       6.00       54.6       27.4       54.3