{Use Pascal Pretty Printer, (C)2001 by Lavrov Dmyry}
unit runge;
interface
type real=double;
const maxlongint=$7fFFffFF;
const MaxRealArraySize=maxlongint div sizeof(real);
type TrealArray=array[1..MaxRealArraySize]of real;pRealArray=^TrealArray;
type Tuserfunc=procedure(t:real;n:longint;v,d:pRealArray);
procedure initRK4(_f:TuserFunc;_n:longint;_dt:real;_p:pointer{pointer to array of real});
procedure StepRK4;
var t:real;
Phase:longint;
dt:real;
implementation
var p:PrealArray;
tmp1,tmp2,tmp3:PrealArray;
f:TuserFunc;
nEQN:longint;

procedure initRK4(_f:TuserFunc;_n:longint;_dt:real;_p:pointer{pointer to array of real});
begin
  (**************disposing*******************)
  if tmp1<>nil then dispose(tmp1);
  if tmp2<>nil then dispose(tmp2);
  if tmp3<>nil then dispose(tmp3);
  nEQN:=_n;
  f:=_f;
  dt:=_dt;
  p:=_p;
  getmem(tmp1,sizeof(real)*nEQN);
  getmem(tmp2,sizeof(real)*nEQN);
  getmem(tmp3,sizeof(real)*nEQN);
end;
procedure StepRK4;
var i:longint;
begin
  if f=nil then exit;
  phase:=0;
  f(T,neqn,p,tmp1); {0}
  for i:=1 to neqn do tmp2^[i]:=p^[i]+dt*tmp1^[i]/2;
  inc(phase);
  f(T+dt/2,neqn,tmp2,tmp3);{1}
  for i:=1 to neqn do begin
    tmp1^[i]:=tmp1^[i]+2*tmp3^[i];
    tmp2^[i]:=p^[i]+dt*tmp3^[i]/2;
  end;
  inc(phase);
  f(T+dt/2,neqn,tmp2,tmp3);{2}
  for i:=1 to neqn do begin
    tmp1^[i]:=tmp1^[i]+2*tmp3^[i];
    tmp2^[i]:=p^[i]+dt*tmp3^[i];
  end;
  inc(phase);
  f(T+dt,neqn,tmp2,tmp3);{3}
  for i:=1 to neqn do begin
   p^[i]:=p^[i]+dt*(tmp1^[i]+tmp3^[i])/6;
  end;
   t:=t+dt;
end;
begin
  tmp1:=nil;
  tmp2:=nil;
  tmp3:=nil;
end.