Изменить программу интегрирования методом второго порядка

Const
Dt=0.1;
V=-0.5;
Q=0.0;
 
Var
f:Text;
V_old,V_new,Q_old,Q_new,t:Real;
i:integer;
 
Begin
Assign(f,'D:\Result1.txt');
Rewrite(f);
t:=0;
Q_old:=Q;
V_old:=V;
Writeln(f,t:5:1,Q:15:3,V:15:3);
for i:=1 to 350 do
begin
t:=Dt*i;
Q_new:=Q_old+V_old*Dt;
V_new:=V_old+Exp(-Q_old)*(Exp(-Q_old)-1)*Dt;
Writeln(f,t:5:1,Q_new:15:3,V_new:15 :3);
Q_old:=Q_new;
V_old:=V_new;
end;
Close(f);
End.

program RungeKutta4;
 
  function f(t, y: real): real;
  begin
    f := exp(-t) * (exp(-t) - 1);
  end;
 
  procedure RungeKutta4Step(var t, y: real; h: real);
  var
    k1, k2, k3, k4: real;
  begin
    k1 := f(t, y);
    k2 := f(t + (h / 2), y + k1 * (h / 2));
    k3 := f(t + (h / 2), y + k2 * (h / 2));
    k4 := f(t + h, y + k3 * h);
    y  := y + (k1 + 2 * k2 + 2 * k3 + k4) * h / 6;
    t  := t + h;
  end;
 
  procedure RungeKutta4(t0, tfin, h, y0: real; var yfin: real; Show: boolean);
  var
    t: real;
    y: real;
  begin
    t := t0;
    y := y0;
    if Show then
      writeln('y(', t: 12: 5, ')=', y: 12: 5);
    while t <= tfin do
    begin
      RungeKutta4Step(t, y, h);
      if Show then
        writeln('y(', t: 12: 5, ')=', y: 12: 5);
    end;
    yfin := y;
  end;
 
  procedure Differ(t0, tfin, y0: real; eps: real);
  var
    h: real;
    y1, y2: real;
  begin
    h := (tfin - t0) / 1;
    RungeKutta4(t0, tfin, h, y0, y2, False);
    repeat
      y1 := y2;
      h  := h / 2;
      RungeKutta4(t0, tfin, h, y0, y2, False);
    until abs(y2 - y1) < eps;
    RungeKutta4(t0, tfin, h, y0, y2, True);
  end;
 
begin
  Differ(0, 35, -0.5, 0.001);
end.

procedure RungeKutta4(t0, tfin, h, y0: real; var yfin: real; Show: boolean);
  var
    t: real;
    y: real;
    f: Text;
  begin
    t := t0;
    y := y0;
    if Show then
    begin
      Assign(f, 'Result.txt');
      Rewrite(f);
    end;
    if Show then
      writeln(f, 'y(', t: 12: 5, ')=', y: 12: 5);
    while t <= tfin do
    begin
      RungeKutta4Step(t, y, h);
      if Show then
        writeln(f, 'y(', t: 12: 5, ')=', y: 12: 5);
    end;
    yfin := y;
    if Show then
      Close(f);
  end;

{
Решение ОДУ 2 порядка модифицированным методом
Эйлера по формуле прогноз+коррекция.
 
Вид ОДУ
y''=f(t, y, y')
 
Обозначения в программе
dy  - y'
ddy - y''
Приводим исходное ОДУ к системе
y1=y
y1'=y2=f1(t, y1, y2)
y2'=y''=f2(t, y1, y2)=f(t, y, y')
 
y1'=f1(t, y1, y2)
y2'=f2(t, y1, y2)
где
f1(t, y1, y2)=y2
f2(t, y1, y2)=f(t, y1, y2)
}
program EulerMethodModificated;
 
  function f(t, y, dy: real): real;
  begin
    f := -sin(t);
  end;
 
  function f1(t, y, dy: real): real;
  begin
    f1 := dy;
  end;
 
  function f2(t, y, dy: real): real;
  begin
    f2 := f(t, y, dy);
  end;
 
  procedure EulerModStep(var t, y1, y2: real; h: real);
  var
    y1p, k1, k2: real;
    y2p, m1, m2: real;
  begin
    y1p := y1 + h * f1(t, y1, y2);
    k1  := f1(t, y1, y2);
    k2  := f1(t, y1p, y2);
    y2p := y2 + h * f2(t, y1, y2);
    m1  := f2(t, y1, y2);
    m2  := f2(t, y1, y2p);
 
    y1 := y1 + (k1 + k2) * h / 2;
    y2 := y2 + (m1 + m2) * h / 2;
    t  := t + h;
  end;
 
  procedure SolveODE2(t0, tfin, h, y0, dy0: real; var yfin: real;
    Show: boolean; Nprint: integer);
  var
    t: real;
    y, dy: real;
    Tprint: real;
  begin
    t  := t0;
    y  := y0;
    dy := dy0;
    if Show then
    begin
      writeln('y(', t: 12: 5, ')=', y: 12: 5);
      Tprint := (tfin - t0) / Nprint + t0;
    end;
    while t <= tfin do
    begin
      EulerModStep(t, y, dy, h);
      if Show and (t >= Tprint) then
      begin
        writeln('y(', t: 12: 5, ')=', y: 12: 5);
        Tprint := Tprint + (tfin - t0) / Nprint;
      end;
    end;
    yfin := y;
  end;
 
  procedure Differ(t0, tfin, y0, dy0: real; eps: real);
  var
    h: real;
    y1, y2: real;
  begin
    h := (tfin - t0) / 10;
    SolveODE2(t0, tfin, h, y0, dy0, y2, False, 0);
    repeat
      y1 := y2;
      h  := h / 2;
      SolveODE2(t0, tfin, h, y0, dy0, y2, False, 0);
    until (abs(y2 - y1) < eps) or (h < eps / 1000);
    SolveODE2(t0, tfin, h, y0, dy0, y2, True, 10);
  end;
 
begin
  Differ(0, 3.5 * 5, 0, 1, 0.01);
end.