#Crib: http://www.math.uni-wuppertal.de/~xsc/literatur/PXSCENGL.pdf #

INT n = 4; # x, xd, y, yd #
[]STRING cmm = ("x ", "xd ", "y ", "yd ");
MODE R = LONG LONG REAL;
MODE RF = PROC(R)R;
RF rsqrt = long long sqrt;
MODE SV = [n]R;
PROC hlf = (SV x)SV: (SV a; FOR j TO n DO a[j] := 0.5*x[j] OD; a);
PROC prs = (STRING m, SV x)VOID: (print((m, newline));
   FOR j TO n DO print((cmm[j], x[j], newline)) OD);
OP * = (R s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a);
OP * = (REAL s, SV v)SV: (SV a; FOR j TO n DO a[j] := s*v[j] OD; a);
OP + = (SV x, y)SV: (SV a; FOR j TO n DO a[j] := x[j] + y[j] OD; a);
# This code is to solve the ordinary differential equation dy/dx = f(x, y) #
# in particular a Newtonian prabolic orbit #
PROC f = (R tau, SV y)SV: (R d2 = y[1]*y[1] + y[3]*y[3], d = rsqrt(d2), dm3 = 1.0/(d*d2);
   (SV a; FOR i BY 2 TO n DO a[i] := y[i+1]; a[i+1] := -y[i]*dm3 OD; a));
PROC step = (R dx, x, SV y)SV: (
  SV k1 = dx*f(x,y);
  SV k2 = dx*f(x+dx/2, y + hlf(k1));
  SV k3 = dx*f(x+dx/2, y + hlf(k2));
  SV k4 = dx*f(x+dx, y+k3);
  y + 1.0/R(6)*(k1+2.0*k2 + 2.0*k3 + k4));

(SV x := (1, 0, 0, rsqrt(2));
PROC rep = (INT j)VOID: (print((j, newline));
    prs("dvs", x); prs("psh", f(0, x));
# This orbit should be 4x + y^2 = 4 #
    print(("parabola ", 4.0*x[1] + x[3]*x[3] - 4.0, newline)));
INT c = 1000, top= 1886;
rep(0);
FOR k TO top DO x := step (1/c, 0, x);
  (k=1 OR k=c OR k=top | rep(k)) OD)

#
         +0
dvs
x +1.0000000000000000000000000000000000000000000000000000000e  +0
xd +0.0000000000000000000000000000000000000000000000000000000e  +0
y +0.0000000000000000000000000000000000000000000000000000000e  +0
yd +1.4142135623730950488016887242096980785696718753769480732e  +0
psh
x +0.0000000000000000000000000000000000000000000000000000000e  +0
xd -1.0000000000000000000000000000000000000000000000000000000e  +0
y +1.4142135623730950488016887242096980785696718753769480732e  +0
yd +0.0000000000000000000000000000000000000000000000000000000e  +0
parabola +0.0000000000000000000000000000000000000000000000000000000e  +0
         +1
dvs
x +9.9999950000016666660416670833330647787898762266032307943e  -1
xd -9.9999933333395833258854253678286075958453244303268312425e  -4
y +1.4142133266709230415893014341028676859988717415379228556e  -3
yd +1.4142128552669031172422413011008233940335901597071605970e  +0
psh
x -9.9999933333395833258854253678286075958453244303268312425e  -4
xd -9.9999800000291666283333812499418653035453832175739217091e  -1
y +1.4142128552669031172422413011008233940335901597071605970e  +0
yd -1.4142112053537614582930891329407642314944199331941946956e  -3
parabola +5.5555555555567925328125014607737670910000000000000000000e -20
      +1000
dvs
x +6.0872178128246237976523436787123492709799790363920578930e  -1
xd -6.3583414768928659846728249324615560946736911510572622206e  -1
y +1.2510447133776198205656091553179136092682256145528563566e  +0
yd +1.0164850878472664560606953966821124000728629618007499753e  +0
psh
x -6.3583414768928659846728249324615560946736911510572622206e  -1
xd -2.2603562033334069797458260336624200471790146961058993350e  -1
y +1.0164850878472664560606953966821124000728629618007499753e  +0
yd -4.6454829866164964822358351853635925166264652897582416296e  -1
parabola -5.9551711237792861314126313566601093844137360000000000000e -14
      +1886
dvs
x -2.7005598369929085610238438330565954847988287501449900003e  -4
xd -7.0710677474212395150648827983061932715799462458924024725e  -1
y +2.0002700377535487478042556518737396274884707942525239701e  +0
yd +7.0701131486847334657277412479760674779817836806990513408e  -1
psh
x -7.0710677474212395150648827983061932715799462458924024725e  -1
xd +3.3743327235504611523351823404266866842181173369466859625e  -5
y +7.0701131486847334657277412479760674779817836806990513408e  -1
yd -2.4993249739820668560100148968070667033844519810679595182e  -1
parabola -2.1383055312957915397238306338359554700526262000000000000e -13
#