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

INT n = 3; # alpha, x, y #
[]STRING cmm = ("alpha ", "x ", "y ");
MODE R = LONG LONG REAL;
MODE RF = PROC(R)R;
RF rsin = long long sin, rcos = long long cos;
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 beam bent so its ends meet #
# 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 f = (R tau, SV y)SV: (-y[3], rcos(y[1]), rsin(y[1]));
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 := (0, 0, 1);
PROC rep = (REAL j)VOID: (print((newline, j, newline));
    prs("dvs", x); (0=1|prs("psh", f(0, x))));
rep(0);
R tau := 0; REAL dt = 0.0001, dr=0.1; REAL tr := 0;
WHILE x[3] > 0 DO x := step (dt, tau, x); tau := tau + dt;
  IF tau > tr-dt/2.0 THEN tr := tr + 0.1; rep(SHORTEN SHORTEN tau) FI OD)