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

INT n = 2;
MODE R = LONG LONG REAL;
MODE RF = PROC(R)R;
RF rsin = long long sin, rcos = long long cos;
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((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 circular harmonic oscillator #
PROC f = (R tau, SV y)SV: (-y[2], 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 := (1, 0);

(INT c = 10000;
prs("ini", x);
TO c DO x := step (1/c, 0, x) OD;
print((c, newline)); prs("end", x);
print((1.0 - (x[1]*x[1] + x[2]*x[2]), 
newline, "cos ", x[1]-rcos(1), newline, "sin ", x[2]-rsin(1))))

#
ini
+1.0000000000000000000000000000000000000000000000000000000e  +0
+0.0000000000000000000000000000000000000000000000000000000e  +0
     +10000
end
+5.4030230586813971810212490461842959905274347886258519194e  -1
+8.4147098480789650620219196619657541787481134949039919280e  -1
+1.3888888871527777777776813368057966579928600000000000000e -22
cos +7.0118829717545299532043305824466296427258410000000000000e -19
sin -4.5031035543372358174775171130797187286931686700000000000e -19
#