comment AP 260;
	  procedure RK5n(x,xa,b,fxj,j,e,d,fi,n,l,pos); value fi,n,l,pos; integer j,n,l;
	  Boolean fi,pos; real b,fxj; array x,xa,e,d;
 	  begin integer i; Boolean first,fir,rej; real fhm,s,s0,cond0,s1,cond1,h,absh,
	    tol,fh,hl,mu,mu1; array y,xl,discr[0:n],k[0:5,0:n],e1[1:2];
	    procedure RKstep(h,d); value h,d; integer d; real h;
	    begin integer i;
	      procedure F(t); value t; integer t;
	      begin integer i; real p; for j:= 0 step 1 until n do y[j]:= fxj; p:=
	        h/sqrt(SUM(i,0,n,y[i]**2)); for i:= 0 step 1 until n do k[t,i]:= y[i] * p
	        end F;
	      if d = 2 then goto integrate; if d = 1 then begin for i:= 0 step 1 
	      until n do k[0,i]:= k[0,i] * mu; goto A end; for i:= 0 step 1 until n
	      do x[i]:= xl[i]; F(0);
	    A:  for i:= 0 step 1 until n do x[i]:= xl[i] + k[0,i]/4.5; F(1); for i:= 0
	      step 1 until n do x[i]:= xl[i] + (k[0,i] + k[1,i] * 3)/12; F(2); for i:= 0
	      step 1 until n do x[i]:= xl[i] + (k[0,i] + k[2,i] * 3)/8; F(3); for i:= 0
	      step 1 until n do x[i]:= xl[i] + (k[0,i] * 53 - k[1,i] * 135 + k[2,i] *
	      126 + k[3,i] * 56)/125; F(4); if d <= 1 then begin for i:= 0 step 1 
	      until n do x[i]:= xl[i] + (k[0,i] * 133 - k[1,i] * 378 + k[2,i] * 276 +
	      k[3,i] * 112 + k[4,i] * 25)/168; F(5); for i:= 0 step 1 until n do
	      discr[i]:= abs(k[0,i] * 21 - k[2,i] * 162 + k[3,i] * 224 - k[4,i] * 125
	      + k[5,i] * 42)/14; goto end end;
	    integrate:  for i:= 0 step 1 until n do x[i]:= xl[i] + (- k[0,i] * 63 + 
	      k[1,i] * 189 - k[2,i] * 36 - k[3,i] * 112 + k[4,i] * 50)/28; F(5); for
	      i:= 0 step 1 until n do x[i]:= xl[i] + (k[0,i] * 35 + k[2,i] * 162 + k[4,i]
	      * 125 + k[5,i] * 14)/336;
	    end:  end RKstep;
	    real procedure fzero;
	    begin if s = s0 then fzero:= cond0 else if s = s1 then fzero:= cond1 
	      else begin RKstep(s-s0,3); fzero:= b end end fzero;
	    if fi then begin for i:= 0 step 1 until n do d[i+3]:= xa[i]; d[1]:= d[2]:= 0
	    end; for i:= 0 step 1 until n do x[i]:= xl[i]:= d[i+3]; s:= d[1]; first:= fir
	    := true; h:= e[0] + e[1]; for i:= 1 step 1 until n do begin absh:= e[2*i]
	    + e[2*i+1]; if h > absh then h:= absh end; if fi then begin j:= l; if fxj
	    * h < 0 == pos then h:= - h end else if d[2] * h < 0 then h:= - h; i:= 0;
	  again:  RKstep(h,i); rej:= false; fhm:= 0; absh:= abs(h); for i:= 0 step 1
	    until n do begin tol:= e[2*i] * abs(k[0,i]) + e[2*i+1] * absh; rej:= tol <
	    discr[i] | rej; fh:= discr[i]/tol; if fh > fhm then fhm:= fh end; mu:= 1/
	    (1 + fhm) + .45; if rej then begin h:= h * mu; i:= 1; goto again end;
	    if first then begin first:= fir; hl:= h; h:= mu * h end else begin fh:=
	    mu * h/hl + mu - mu1; hl:= h; h:= fh * h end;
	  accept:  RKstep(hl,2); mu1:= mu; s:= s + hl; if fir then begin cond0:= b;
	    fir:= false; if ! fi then h:= d[2] end else begin d[2]:= h; cond1:= b; if
	    cond0 * cond1 <= 0 then goto zero; cond0:= cond1 end; for i:= 0 step 1
	    until n do d[i+3]:= xl[i]:= x[i]; d[1]:= s0:= s; i:= 0; goto again;
	  zero:  e1[1]:= e[2*n+2]; e1[2]:= e[2*n+3]; s1:= s; s:= ZERO(s,s0,s1,fzero,
	    e1); RKstep(s-s0,3); for i:= 0 step 1 until n do d[i+3]:= x[i]; d[1]:= s
	    end RK5n;