then k[0,i]:= p * y[i] end end
	      else for i:= 0 step 1 until n do begin if i != iv then k[0,i]:= k[0,i] *
	      mu end; for i:= 0 step 1 until n do x[i]:= xl[i] + (if i = iv then h
	      else k[0,i])/4.5; F(1); for i:= 0 step 1 until n do x[i]:= xl[i] + (if
	      i = iv then h * 4 else (k[0,i] + k[1,i] * 3))/12; F(2); for i:= 0 step
	      1 until n do x[i]:= xl[i] + (if i = iv then h * .5 else (k[0,i] + k[2,i]
	      * 3)/8); F(3); for i:= 0 step 1 until n do x[i]:= xl[i] + (if i = iv
	      then h * .8 else (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] + (if i = iv then h else (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 begin if i != iv then discr[i]:= abs(k[0,i] * 21 - k[2,i] *
	      162 + k[3,i] * 224 - k[4,i] * 125 + k[5,i] * 42)/14 end; goto end end;
	    integrate:  for i:= 0 step 1 until n do x[i]:= xl[i] + (if i = iv then h
	      else (- 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 begin if i != iv then x[i]
	      := xl[i] + (k[0,i] * 35 + k[2,i] * 162 + k[4,i] * 125 + k[5,i] * 14)/336
	      end;
	    end:  end RKstep;
	    real procedure fzero;
	    begin if s = x0 then fzero:= cond0 else if s = x1 then fzero:= cond1
	      else begin RKstep(s-xl[iv],3); fzero:= b end end fzero;
	    if fi then begin for i:= 0 step 1 until n do d[i+3]:= xa[i]; d[0]:= d[2]:=
	    0 end; d[1]:= 0; for i:= 0 step 1 until n do x[i]:= xl[i]:= d[i+3]; iv:=
	    d[0]; h:= d[2]; first:= fir:= true; y[0]:= 1; goto change;
	  again:  absh:= abs(h); if absh < hmin then begin h:= if h > 0 then hmin
	    else - hmin; absh:= abs(h) end; RKstep(h,i); rej:= false; fhm:= 0;
	    for i:= 0 step 1 until n do begin if i!=iv then 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 end; mu:= 1/(1 + fhm) + .45; if rej
	    then begin if absh <= hmin then begin for i:= 0 step 1 until n do begin
	    if i!=iv then x[i]:= xl[i] + k[0,i] else x[i]:= xl[i] + h end; d[1]:= d[1] + 1;
	    first:= true; goto next end; h:= h * mu; i:= 0; goto again end; if first
	    then begin first:= fir; hl:= h; h:= mu * h; goto accept end; fh:= mu *
	    h/hl + mu - mu1; hl:= h; h:= fh * h;
	  accept:  RKstep(hl,2); mu1:= mu;
	  next:  if fir then begin fir:= false; cond0:= b; if ! (fi | rej) 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];
	  change:  iv0:= iv; for j:= 1 step 1 until n do y[j]:= fxj; max:= abs(y[iv]);
	    for i:= 0 step 1 until n do begin if abs(y[i]) > max then begin max:=
	    abs(y[i]); iv:= i end end; if iv0!=iv then begin first:= true; d[0]:= iv;
	    d[2]:= h:= y[iv]/y[iv0] * h end; x0:= xl[iv]; if fir then begin hmin:=
	    e[0] + e[1]; for i:= 1 step 1 until n do begin h:= e[2*i] + e[2*i+1]; if
	    h < hmin then hmin:= h end; h:= e[2*iv] + e[2*iv+1]; if (fi & (y[l]/y[iv]
	    * h < 0 == pos)) | (! fi & d[2] * h < 0) then h:= - h end; i:= 1; goto
	    again;
	  zero: e1[1]:= e[2*n+2]; e1[2]:= e[2*n+3]; x1:= x[iv]; x0:= ZERO(s,x0,x1,
	    fzero,e1); RKstep(x0-xl[iv],3); for i:= 0 step 1 until n do d[i+3]:= x[i]
	    end RK4n;]