Решение систем дифференциальных уравнений

program project;
const jm=5;
type mas=array[1..5]of real;
var J1,J2,C1,C2,Kc,M,J,K1,K2:real;
t0,t,tk,h:real;
q,n:integer;
f,x0,x:mas;
d:text;
l:string;
procedure rp(t:real;var x0,f:mas);
begin
f[1]:=(1/J1)*(x[3]-(x[1]-x[2])*Kc); {DWj1/Dt}
f[2]:=(1/J2)*((x[1]-x[2])*Kc-x[4]); {DWj2/Dt}
f[3]:=C1*((M-x[3])/K1 - x[1]); {DMc1/Dt}
f[4]:=C2*(x[2]-x[5]); {DMc2/Dt}
f[5]:=1/J*(x[4]-x[5]/K2); {DJ/Dt}
end;
procedure eu(var t,t0,h:real;var x0:mas;var x:mas);
var j:integer;
begin
rp(t,x0,f);
for j:=1 to 5 do
x[j]:=x0[j]+h*f[j];
t:=t0+h;
end;
begin
writeln('ishodnie dannie');
readln(l);
assign(d,l);
reset(d);
read(d,J1,J2,C1,C2,Kc,M,J,K1,K2);
close(d) ;
writeln('rezultati');
readln(l);
assign(d,l);
rewrite(d);
t0:=0.0; tk:=5.0;
h:=0.01;
n:=0;
x0[1]:=0.0; x0[2]:=0.0; x0[3]:=0.0; x0[4]:=0.0; x0[5]:=0.0;
t:=t0;
writeln('t=',t:2:2,'; x1=',x0[1]:5:6,'; x2=',x0[2]:5:6,'; x3=',x0[3]:5:6,'; x4=',x0[4]:5:6,'; x5=',x0[5]:5:6,';');
writeln(d,'t=',t:2:2,'; x1=',x0[1]:5:6,'; x2=',x0[2]:5:6,'; x3=',x0[3]:5:6,'; x4=',x0[4]:5:6,'; x5=',x0[5]:5:6,';');
while t<=tk do
begin
eu(t,t0,h,x0,x);
n:=n+1;
if n=10 then
begin
writeln(d,'t=',t:2:2,'; x1=',x[1]:5:6,'; x2=',x[2]:5:6,'; x3=',x[3]:5:6,'; x4=',x[4]:5:6,'; x5=',x[5]:5:6,';');
writeln('t=',t:2:2,'; x1=',x[1]:5:6,'; x2=',x[2]:5:6,'; x3=',x[3]:5:6,'; x4=',x[4]:5:6,'; x5=',x[5]:5:6,';');
n:=0;
end;
for q:=1 to 5 do
x0[q]:=x[q];
t0:=t;
end;
close(d);
end.