#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 = 10000, top= 55000; R dx = 1/c;
rep(0);
FOR k TO top DO x := step (dx, 0, x);
  (k=1 OR (k MOD c) = 0 OR k=top | rep(k)) OD)

#         +0
dvs
x +1.00000000000000000000000000000000000000000000000000000000000000e  +0
xd +0.00000000000000000000000000000000000000000000000000000000000000e  +0
y +0.00000000000000000000000000000000000000000000000000000000000000e  +0
yd +1.41421356237309504880168872420969807856967187537694807317667974e  +0
psh
x +0.00000000000000000000000000000000000000000000000000000000000000e  +0
xd -1.00000000000000000000000000000000000000000000000000000000000000e  +0
y +1.41421356237309504880168872420969807856967187537694807317667974e  +0
yd +0.00000000000000000000000000000000000000000000000000000000000000e  +0
parabola +0.00000000000000000000000000000000000000000000000000000000000000e  +0
         +1
dvs
x +9.99999995000000016666666604166667083333330647786472737630132853e  -1
xd -9.99999993333333395833332588541675367838442325847475158678922501e  -5
y +1.41421356001607245368536503017788742465402489787791994205820225e  -4
yd +1.41421355530202729586177791616761541682110383262411492058686998e  +0
psh
x -9.99999993333333395833332588541675367838442325847475158678922501e  -5
xd -9.99999980000000291666662833333381249999418652350667046359779982e  -1
y +1.41421355530202729586177791616761541682110383262411492058686998e  +0
yd -1.41421353880286933628698666311526347883250219011089942454470199e  -4
parabola +5.55555555555555567925347031250001460780000000000000000000000000e -26
     +10000
dvs
x +6.08721781282468751692660505803174200265888275298245741252329202e  -1
xd -6.35834147689268604754633457052697722842831821201157217756724281e  -1
y +1.25104471337763343240798572367813295315942208423974820689400843e  +0
yd +1.01648508784727860505754688668825036381473443944872199550603365e  +0
psh
x -6.35834147689268604754633457052697722842831821201157217756724281e  -1
xd -2.26035620333335739558702594732548478187805279378849917518198682e  -1
y +1.01648508784727860505754688668825036381473443944872199550603365e  +0
yd -4.64548298661639649391653880805302416042850401009555680484639839e  -1
parabola -5.95464417034105218716849723792691633217301774400000000000000000e -18
     +20000
dvs
x -8.08594603928763757933625488764077898515421823920041196779860918e  -2
xd -7.06572714825347830630488227851612025471746262620298675203375059e  -1
y +2.07928782076255750847655386436364664839057631788351550655651463e  +0
yd +6.79629542163354297188820805010140567709913083027041187376821185e  -1
psh
x -7.06572714825347830630488227851612025471746262620298675203375059e  -1
xd +8.97434148374444182264227070510602640848572109056195159360225951e  -3
y +6.79629542163354297188820805010140567709913083027041187376821185e  -1
yd -2.30773725867677725146277617755425765513526988514796059616794295e  -1
parabola -2.37422658302811365871405444748382674333110275880000000000000000e -17
     +30000
dvs
x -7.75726623466793165272655302622630841738590072468473895073458386e  -1
xd -6.78932126976413522689937460140925159651921588943738147538280182e  -1
y +2.66512785694554860949873044103396424209322859836076105579762108e  +0
yd +5.09493100083029014757365854654172017849579191784035050717589420e  -1
psh
x -6.78932126976413522689937460140925159651921588943738147538280182e  -1
xd +3.62725947696178928252981859682932813239334441397156109197553253e  -2
y +5.09493100083029014757365854654172017849579191784035050717589420e  -1
yd -1.24620065677485785192611698602253548938331536247330496909699995e  -1
parabola -4.74422936474959736487132028423917795430874115820000000000000000e -17
     +40000
dvs
x -1.43628491597273211684678882829342515976918628469014727463476569e  +0
xd -6.42376833488174989335804721722704890545573918338482975748293954e  -1
y +3.12172062553504752155487112627847502426474078163447188499339863e  +0
yd +4.11553057140130699773201921208241375985543944699531326779816550e  -1
psh
x -6.42376833488174989335804721722704890545573918338482975748293954e  -1
xd +3.53975417216011261571747975754473095748165772150660257882349804e  -2
y +4.11553057140130699773201921208241375985543944699531326779816550e  -1
yd -7.69354567862477486057116796756975410235194880031846863823745653e  -2
parabola -7.53149329346107004030000776205826127273329123430000000000000000e -17
     +50000
dvs
x -2.06170354394960126602381594313020227691951010114206100324377622e  +0
xd -6.09239908725110659833712053086914023276163104929951802120670241e  -1
y +3.49954485266275805132799629698339204891903173440772215224831487e  +0
yd +3.48182369065250270162744478431629570663697716504174329694254389e  -1
psh
x -6.09239908725110659833712053086914023276163104929951802120670241e  -1
xd +3.07681618748136252044925044702652318158812779206652061878570544e  -2
y +3.48182369065250270162744478431629570663697716504174329694254389e  -1
yd -5.22260161170536662975510253291276874320869652500730161475779558e  -2
parabola -1.06363133128754140584937939304370127584785134320000000000000000e -16
     +55000
dvs
x -2.36257531466267626833087516671515017160637782612458094710346185e  +0
xd -5.94440061483565566965427277705964364264949417541064577397498494e  -1
y +3.66746523618843658568955312785486728890940354584215292865142018e  +0
yd +3.24169432128748264472828432904265944027938932771592227781244397e  -1
psh
x -5.94440061483565566965427277705964364264949417541064577397498494e  -1
xd +2.84548859403596451176808658424462793898161407594436910605428590e  -2
y +3.24169432128748264472828432904265944027938932771592227781244397e  -1
yd -4.41709960898647730905265350913433630999084398489011679767589481e  -2
parabola -1.22874049950468686316834336937520843658627949764000000000000000e -16
#