Двухстадийный метод Розенброка третьего порядка точности

Метод применяется для решения систем уравнений, в которых правая часть не зависит от t, т.е. для систем вида:

(1)        dy/dt = f (y);   y ( t₀ ) = y₀

Ограничение на вид системы не является существенным, для систем вида

dy/dt = f (t, y);   y ( t₀) = y₀

всегда можно ввести еще одну переменную y_n+1 = t и добавить еще одно уравнение

dy_n+1/dt = 1;

тогда, заменив в правой части вхождения t на y_n+1, получим систему вида (1).

---

Схема метода

y_n+1 = y_n + w₁k₁ + w₂k₂

k₁ = h[I - ha₁ J(y_n )]^-1f(y_n )

k₂ = h[I - ha₂ J(y_n +c₁k₁)]^-1f(y_n + b₁k₁)


J - якобиан, I - единичная матрица.

---

Параметры схемы:

a₁ = 1 + ||/ 6 /6 = 1.40824829046386;
a₂ = 1 - ||/ 6 /6 = 5.91751709536137e-1;
w₁ = -0.41315432;
w₂ = 1.41315432;
b₁ = c₁ = ( - 6 - ||/ 6  + [58 + 20||/ 6  ]^1/2 )/( 6 + 2||/ 6  ) =
= 1.73786673924946e-1;

---

Метод у нас записан немножко по другому, мы находим k₁ и k₂ решая системы уравнений

[I - ha₁ J(y_n )] k₁ = hf ( y_n )

[I - ha₂ J(y_n +c₁k₁)]k₂ = h f ( y_n + b₁k₁)

---

Текст программы на языке Паскаль. Программа несколько корявая, но она работает, я облазил много сайтов, везде или метод криво изложен или в тексте программ ошибки(?), по крайней мере у меня эти программы не работали. Программа писалась под те задачи, которые у нас решаются, поэтому для произвольной системы уравнений может понадобиться некоторая переделка. Например, может понадобиться переписать процедуру jacob, у нас она работает только для положительных y[i] (отрицательные у нас невозможны по условиям задачи). Программа только для системы из пяти уравнений, но необходимые изменения (при желании) внести будет нетрудно. Что касается выбора шага, то об этом тоже надо позаботиться самому.

{Semi-Implicit forms of Rosenbrock
These are suitable for autonomous system of equations
of the form
dy/dt = f(y), y(t0)=y0 }

{$n+}{directiva compilyatora dlya vychisleniy s 
dvoinoi tochnostyu}
TYPE
   VTAB = array[1..5] of double;
   ATAB = array[1..5,1..5] of double;
var
y,ynext:vtab;
a,jcb:atab;
eps,eps1,h,x,xnext:double;
buffer:double;
i,j,k:integer;
res:text;{TEXTOVYI FAIL, V KOTORYJ VYVODITSYA RESULTAT}

{--------------------}
procedure dy(var y,f:vtab);
{PROCEDURA DLYA NAHOZHDENIYA dy/dt
f[i] - proizvodnye dy[i]/dt,
SJUDA VCTAVITE SVOI URAVNENIYA

   SLEDI ZA PEREDACHEI ZNACHENIJ ! !

}
begin
f[1]:=y[1];
f[2]:=y[2];
f[3]:=y[3];
f[4]:=y[4];
f[5]:=y[5];
end;
{VNIMANIE! V PROTSEDURE JACOB Y[i] DOLZHNY BYT' POLOZHITELNY}
procedure jacob(var y:vtab;var jcb:atab);
{VNIMANIE! V PROTSEDURE JACOB Y[i] DOLZHNY BYT' POLOZHITELNY;
PRI NEOBHODIMOSTI SAMI POPRAVITE}
var x0,x1,delta:double;
    y0,y1,f0,f1:vtab;
    k:integer;
begin
delta:=1e-7;

for k:=1 to 5 do begin
y0[k]:=y[k];
y1[k]:=y[k];
end;
y0[1]:=y0[1]*(1-delta);
y1[1]:=y1[1]*(1+delta);
dy(y0,f0);
dy(y1,f1);
jcb[1,1]:=(f1[1]-f0[1])/(2*delta*y[1]);{d_f1 po d_y1}
jcb[2,1]:=(f1[2]-f0[2])/(2*delta*y[1]);{d_f2 po d_y1}
jcb[3,1]:=(f1[3]-f0[3])/(2*delta*y[1]);{d_f3 po d_y1}
jcb[4,1]:=(f1[4]-f0[4])/(2*delta*y[1]);{d_f4 po d_y1}
jcb[5,1]:=(f1[5]-f0[5])/(2*delta*y[1]);{d_f5 po d_y1}

y0[1]:=y[1];y1[1]:=y[1];{vozvraschaem znacheniya}
y0[2]:=y0[2]*(1-delta);
y1[2]:=y1[2]*(1+delta);
dy(y0,f0);
dy(y1,f1);
jcb[1,2]:=(f1[1]-f0[1])/(2*delta*y[2]); {d_f1 po d_y2}
jcb[2,2]:=(f1[2]-f0[2])/(2*delta*y[2]); {d_f2 po d_y2}
jcb[3,2]:=(f1[3]-f0[3])/(2*delta*y[2]); {d_f3 po d_y2}
jcb[4,2]:=(f1[4]-f0[4])/(2*delta*y[2]); {d_f4 po d_y2}
jcb[5,2]:=(f1[5]-f0[5])/(2*delta*y[2]); {d_f5 po d_y2}

y0[2]:=y[2];y1[2]:=y[2];{vozvraschaem znacheniya}
y0[3]:=y0[3]*(1-delta);
y1[3]:=y1[3]*(1+delta);
dy(y0,f0);
dy(y1,f1);
jcb[1,3]:=(f1[1]-f0[1])/(2*delta*y[3]); {d_f1 po d_y3}
jcb[2,3]:=(f1[2]-f0[2])/(2*delta*y[3]); {d_f2 po d_y3}
jcb[3,3]:=(f1[3]-f0[3])/(2*delta*y[3]); {d_f3 po d_y3}
jcb[4,3]:=(f1[4]-f0[4])/(2*delta*y[3]); {d_f4 po d_y3}
jcb[5,3]:=(f1[5]-f0[5])/(2*delta*y[3]); {d_f5 po d_y3}

y0[3]:=y[3];y1[3]:=y[3];{vozvraschaem znacheniya}
y0[4]:=y0[4]*(1-delta);
y1[4]:=y1[4]*(1+delta);
dy(y0,f0);
dy(y1,f1);
jcb[1,4]:=(f1[1]-f0[1])/(2*delta*y[4]); {d_f1 po d_y4}
jcb[2,4]:=(f1[2]-f0[2])/(2*delta*y[4]); {d_f2 po d_y4}
jcb[3,4]:=(f1[3]-f0[3])/(2*delta*y[4]); {d_f3 po d_y4}
jcb[4,4]:=(f1[4]-f0[4])/(2*delta*y[4]); {d_f4 po d_y4}
jcb[5,4]:=(f1[5]-f0[5])/(2*delta*y[4]); {d_f5 po d_y4}

y0[4]:=y[4];y1[4]:=y[4];{vozvraschaem znacheniya}
y0[5]:=y0[5]*(1-delta);
y1[5]:=y1[5]*(1+delta);
dy(y0,f0);
dy(y1,f1);
jcb[1,5]:=(f1[1]-f0[1])/(2*delta*y[5]); {d_f1 po d_y5}
jcb[2,5]:=(f1[2]-f0[2])/(2*delta*y[5]); {d_f2 po d_y5}
jcb[3,5]:=(f1[3]-f0[3])/(2*delta*y[5]); {d_f3 po d_y5}
jcb[4,5]:=(f1[4]-f0[4])/(2*delta*y[5]); {d_f4 po d_y5}
jcb[5,5]:=(f1[5]-f0[5])/(2*delta*y[5]); {d_f5 po d_y5}


end;
procedure mult_a_v(var a:atab;var v,c:vtab);
{protsedura umnozhaet matritsu A na vektor v
rezultat  - vektor c}
begin
for i:=1 to 5 do
c[i]:=a[i,1]*v[1]+a[i,2]*v[2]+a[i,3]*v[3]+
      a[i,4]*v[4]+a[i,5]*v[5];
end;
{-------------------}
function det5(var a:atab):double;
{FUNKCIYA VYCHISLYAET OPREDELITEL MATRITSY 5-go PORYADKA a}
begin
det5:=
+a[1,1]*a[2,2]*a[3,3]*a[4,4]*a[5,5]-a[1,2]*a[2,1]*a[3,3]*a[4,4]*a[5,5]
+a[1,3]*a[2,1]*a[3,2]*a[4,4]*a[5,5]-a[1,1]*a[2,3]*a[3,2]*a[4,4]*a[5,5]
-a[1,3]*a[2,2]*a[3,1]*a[4,4]*a[5,5]+a[1,2]*a[2,3]*a[3,1]*a[4,4]*a[5,5]
-a[1,4]*a[2,1]*a[3,2]*a[4,3]*a[5,5]+a[1,1]*a[2,4]*a[3,2]*a[4,3]*a[5,5]
-a[1,1]*a[2,2]*a[3,4]*a[4,3]*a[5,5]+a[1,4]*a[2,2]*a[3,1]*a[4,3]*a[5,5]
-a[1,2]*a[2,4]*a[3,1]*a[4,3]*a[5,5]+a[1,2]*a[2,1]*a[3,4]*a[4,3]*a[5,5]
-a[1,4]*a[2,3]*a[3,1]*a[4,2]*a[5,5]+a[1,3]*a[2,4]*a[3,1]*a[4,2]*a[5,5]
-a[1,3]*a[2,1]*a[3,4]*a[4,2]*a[5,5]+a[1,4]*a[2,1]*a[3,3]*a[4,2]*a[5,5]
-a[1,1]*a[2,4]*a[3,3]*a[4,2]*a[5,5]+a[1,1]*a[2,3]*a[3,4]*a[4,2]*a[5,5]
+a[1,4]*a[2,3]*a[3,2]*a[4,1]*a[5,5]-a[1,3]*a[2,4]*a[3,2]*a[4,1]*a[5,5]
+a[1,3]*a[2,2]*a[3,4]*a[4,1]*a[5,5]-a[1,4]*a[2,2]*a[3,3]*a[4,1]*a[5,5]
+a[1,2]*a[2,4]*a[3,3]*a[4,1]*a[5,5]-a[1,2]*a[2,3]*a[3,4]*a[4,1]*a[5,5]
+a[1,5]*a[2,1]*a[3,2]*a[4,3]*a[5,4]-a[1,1]*a[2,5]*a[3,2]*a[4,3]*a[5,4]
+a[1,1]*a[2,2]*a[3,5]*a[4,3]*a[5,4]-a[1,1]*a[2,2]*a[3,3]*a[4,5]*a[5,4]
-a[1,5]*a[2,2]*a[3,1]*a[4,3]*a[5,4]+a[1,2]*a[2,5]*a[3,1]*a[4,3]*a[5,4]
-a[1,2]*a[2,1]*a[3,5]*a[4,3]*a[5,4]+a[1,2]*a[2,1]*a[3,3]*a[4,5]*a[5,4]
+a[1,5]*a[2,3]*a[3,1]*a[4,2]*a[5,4]-a[1,3]*a[2,5]*a[3,1]*a[4,2]*a[5,4]
+a[1,3]*a[2,1]*a[3,5]*a[4,2]*a[5,4]-a[1,3]*a[2,1]*a[3,2]*a[4,5]*a[5,4]
-a[1,5]*a[2,1]*a[3,3]*a[4,2]*a[5,4]+a[1,1]*a[2,5]*a[3,3]*a[4,2]*a[5,4]
-a[1,1]*a[2,3]*a[3,5]*a[4,2]*a[5,4]+a[1,1]*a[2,3]*a[3,2]*a[4,5]*a[5,4]
-a[1,5]*a[2,3]*a[3,2]*a[4,1]*a[5,4]+a[1,3]*a[2,5]*a[3,2]*a[4,1]*a[5,4]
-a[1,3]*a[2,2]*a[3,5]*a[4,1]*a[5,4]+a[1,3]*a[2,2]*a[3,1]*a[4,5]*a[5,4]
+a[1,5]*a[2,2]*a[3,3]*a[4,1]*a[5,4]-a[1,2]*a[2,5]*a[3,3]*a[4,1]*a[5,4]
+a[1,2]*a[2,3]*a[3,5]*a[4,1]*a[5,4]-a[1,2]*a[2,3]*a[3,1]*a[4,5]*a[5,4]
-a[1,5]*a[2,4]*a[3,1]*a[4,2]*a[5,3]+a[1,4]*a[2,5]*a[3,1]*a[4,2]*a[5,3]
-a[1,4]*a[2,1]*a[3,5]*a[4,2]*a[5,3]+a[1,4]*a[2,1]*a[3,2]*a[4,5]*a[5,3]
+a[1,5]*a[2,1]*a[3,4]*a[4,2]*a[5,3]-a[1,1]*a[2,5]*a[3,4]*a[4,2]*a[5,3]
+a[1,1]*a[2,4]*a[3,5]*a[4,2]*a[5,3]-a[1,1]*a[2,4]*a[3,2]*a[4,5]*a[5,3]
-a[1,5]*a[2,1]*a[3,2]*a[4,4]*a[5,3]+a[1,1]*a[2,5]*a[3,2]*a[4,4]*a[5,3]
-a[1,1]*a[2,2]*a[3,5]*a[4,4]*a[5,3]+a[1,1]*a[2,2]*a[3,4]*a[4,5]*a[5,3]
+a[1,5]*a[2,4]*a[3,2]*a[4,1]*a[5,3]-a[1,4]*a[2,5]*a[3,2]*a[4,1]*a[5,3]
+a[1,4]*a[2,2]*a[3,5]*a[4,1]*a[5,3]-a[1,4]*a[2,2]*a[3,1]*a[4,5]*a[5,3]
-a[1,5]*a[2,2]*a[3,4]*a[4,1]*a[5,3]+a[1,2]*a[2,5]*a[3,4]*a[4,1]*a[5,3]
-a[1,2]*a[2,4]*a[3,5]*a[4,1]*a[5,3]+a[1,2]*a[2,4]*a[3,1]*a[4,5]*a[5,3]
+a[1,5]*a[2,2]*a[3,1]*a[4,4]*a[5,3]-a[1,2]*a[2,5]*a[3,1]*a[4,4]*a[5,3]
+a[1,2]*a[2,1]*a[3,5]*a[4,4]*a[5,3]-a[1,2]*a[2,1]*a[3,4]*a[4,5]*a[5,3]
-a[1,5]*a[2,4]*a[3,3]*a[4,1]*a[5,2]+a[1,4]*a[2,5]*a[3,3]*a[4,1]*a[5,2]
-a[1,4]*a[2,3]*a[3,5]*a[4,1]*a[5,2]+a[1,4]*a[2,3]*a[3,1]*a[4,5]*a[5,2]
+a[1,5]*a[2,3]*a[3,4]*a[4,1]*a[5,2]-a[1,3]*a[2,5]*a[3,4]*a[4,1]*a[5,2]
+a[1,3]*a[2,4]*a[3,5]*a[4,1]*a[5,2]-a[1,3]*a[2,4]*a[3,1]*a[4,5]*a[5,2]
-a[1,5]*a[2,3]*a[3,1]*a[4,4]*a[5,2]+a[1,3]*a[2,5]*a[3,1]*a[4,4]*a[5,2]
-a[1,3]*a[2,1]*a[3,5]*a[4,4]*a[5,2]+a[1,3]*a[2,1]*a[3,4]*a[4,5]*a[5,2]
+a[1,5]*a[2,4]*a[3,1]*a[4,3]*a[5,2]-a[1,4]*a[2,5]*a[3,1]*a[4,3]*a[5,2]
+a[1,4]*a[2,1]*a[3,5]*a[4,3]*a[5,2]-a[1,4]*a[2,1]*a[3,3]*a[4,5]*a[5,2]
-a[1,5]*a[2,1]*a[3,4]*a[4,3]*a[5,2]+a[1,1]*a[2,5]*a[3,4]*a[4,3]*a[5,2]
-a[1,1]*a[2,4]*a[3,5]*a[4,3]*a[5,2]+a[1,1]*a[2,4]*a[3,3]*a[4,5]*a[5,2]
+a[1,5]*a[2,1]*a[3,3]*a[4,4]*a[5,2]-a[1,1]*a[2,5]*a[3,3]*a[4,4]*a[5,2]
+a[1,1]*a[2,3]*a[3,5]*a[4,4]*a[5,2]-a[1,1]*a[2,3]*a[3,4]*a[4,5]*a[5,2]
+a[1,5]*a[2,4]*a[3,3]*a[4,2]*a[5,1]-a[1,4]*a[2,5]*a[3,3]*a[4,2]*a[5,1]
+a[1,4]*a[2,3]*a[3,5]*a[4,2]*a[5,1]-a[1,4]*a[2,3]*a[3,2]*a[4,5]*a[5,1]
-a[1,5]*a[2,3]*a[3,4]*a[4,2]*a[5,1]+a[1,3]*a[2,5]*a[3,4]*a[4,2]*a[5,1]
-a[1,3]*a[2,4]*a[3,5]*a[4,2]*a[5,1]+a[1,3]*a[2,4]*a[3,2]*a[4,5]*a[5,1]
+a[1,5]*a[2,3]*a[3,2]*a[4,4]*a[5,1]-a[1,3]*a[2,5]*a[3,2]*a[4,4]*a[5,1]
+a[1,3]*a[2,2]*a[3,5]*a[4,4]*a[5,1]-a[1,3]*a[2,2]*a[3,4]*a[4,5]*a[5,1]
-a[1,5]*a[2,4]*a[3,2]*a[4,3]*a[5,1]+a[1,4]*a[2,5]*a[3,2]*a[4,3]*a[5,1]
-a[1,4]*a[2,2]*a[3,5]*a[4,3]*a[5,1]+a[1,4]*a[2,2]*a[3,3]*a[4,5]*a[5,1]
+a[1,5]*a[2,2]*a[3,4]*a[4,3]*a[5,1]-a[1,2]*a[2,5]*a[3,4]*a[4,3]*a[5,1]
+a[1,2]*a[2,4]*a[3,5]*a[4,3]*a[5,1]-a[1,2]*a[2,4]*a[3,3]*a[4,5]*a[5,1]
-a[1,5]*a[2,2]*a[3,3]*a[4,4]*a[5,1]+a[1,2]*a[2,5]*a[3,3]*a[4,4]*a[5,1]
-a[1,2]*a[2,3]*a[3,5]*a[4,4]*a[5,1]+a[1,2]*a[2,3]*a[3,4]*a[4,5]*a[5,1];
end;
{-------------------}

procedure kramer(var sys:atab;var d,reshenie:vtab);
var y1,y2,y3,y4,y5,detsys:extended;sys1:atab;
    i4,j4:integer;
{      d - stolbets pravyh chlenov
     sys - matritsa sistemy
reshenie - vektor reshenie}
begin
{ Zapominaem matritsu sistemy }
for i4:=1 to 5 do
for J4:=1 to 5 do
sys1[i4,j4]:=sys[i4,j4];
{--- nahodim opredelitel' matritsy sistemy -------}
detsys:=det5(sys);
if detsys=0 then begin
       writeln(' DETERMINANT = ',detsys);
       writeln;
       writeln('   PRESS   ');readln;halt;
       end;
{---   nahodim reshenie[1]  --------}
   { Zapominaem matritsu sistemy }
        for i4:=1 to 5 do
        for j4:=1 to 5 do
        sys1[i4,j4]:=sys[i4,j4];
   {podstavlyaem vmesto 1-go stolbtsa
        stolbets svobodnyh chlenov}
        for i4:=1 to 5 do
        sys1[i4,1]:=d[i4];
   {nahodim opredelitel' }
       y1:=det5(sys1);
   {nahodim reshenie[1]  }
       reshenie[1]:=y1/detsys;
{---   nahodim reshenie[2]  --------}
   { Zapominaem matritsu sistemy }
        for i4:=1 to 5 do
        for j4:=1 to 5 do
        sys1[i4,j4]:=sys[i4,j4];
   {podstavlyaem vmesto 2-go stolbtsa
        stolbets svobodnyh chlenov}
        for i4:=1 to 5 do
        sys1[i4,2]:=d[i4];
   {nahodim opredelitel' }
       y2:=det5(sys1);
   {nahodim reshenie[2]  }
       reshenie[2]:=y2/detsys;
{---   nahodim reshenie[3]  --------}
   { Zapominaem matritsu sistemy }
        for i4:=1 to 5 do
        for j4:=1 to 5 do
        sys1[i4,j4]:=sys[i4,j4];
   {podstavlyaem vmesto 3-go stolbtsa
        stolbets svobodnyh chlenov}
        for i4:=1 to 5 do
        sys1[i4,3]:=d[i4];
   {nahodim opredelitel' }
       y3:=det5(sys1);
   {nahodim reshenie[3]  }
       reshenie[3]:=y3/detsys;
{---   nahodim reshenie[4]  --------}
   { Zapominaem matritsu sistemy }
        for i4:=1 to 5 do
        for j4:=1 to 5 do
        sys1[i4,j4]:=sys[i4,j4];
   {podstavlyaem vmesto 4-go stolbtsa
        stolbets svobodnyh chlenov}
        for i4:=1 to 5 do
        sys1[i4,4]:=d[i4];
   {nahodim opredelitel' }
       y4:=det5(sys1);
   {nahodim reshenie[4]  }
       reshenie[4]:=y4/detsys;
{---   nahodim reshenie[5]  --------}
   { Zapominaem matritsu sistemy }
        for i4:=1 to 5 do
        for j4:=1 to 5 do
        sys1[i4,j4]:=sys[i4,j4];
   {podstavlyaem vmesto 5-go stolbtsa
        stolbets svobodnyh chlenov}
        for i4:=1 to 5 do
        sys1[i4,5]:=d[i4];
   {nahodim opredelitel' }
       y5:=det5(sys1);
   {nahodim reshenie[5]  }
       reshenie[5]:=y5/detsys;

end;

{-------------------}

procedure rsnbrk(var h,x,xnext:double;var y,ynext:vtab;var jcb:atab);
CONST
   w1=-0.41315432;
   w2=1.41315432;
   a1=1.40824829046386;
   a2=5.91751709536137e-1;
   b1=1.73786673924946e-1;
   c1=1.73786673924946e-1;
VAR y1,k1,k2,f1:vtab;
    j1,j2:atab;
    ir,jr:integer;
begin
jacob(y,jcb);

{POLUCHAEM MATRITSU SISTEMY DLYA NAHOZHDENIYA k1}
for ir:=1 to 5 do
begin
for jr:=1 to 5 do
j1[ir,jr]:=-h*a1*jcb[ir,jr];
j1[ir,ir]:=1-h*a1*jcb[ir,ir];
end;



{NAHODIM f(y_n)}
dy(y,f1);
{POLUCHAEM PRAVUYU CHAST'}
for ir:=1 to 5 do
f1[ir]:=h*f1[ir];
{RESHAEM SISTEMY}
kramer(j1,f1,k1);
{ --------- k1 NAIDENO ----------- }

{ --------- NAHODIM k2 ------------}
for ir:=1 to 5 do
y1[ir]:=y[ir]+c1*k1[ir];
jacob(y1,jcb);
{POLUCHAEM MATRITSU SISTEMY DLYA NAHOZHDENIYA k2}
for ir:=1 to 5 do
begin
for jr:=1 to 5 do
j1[ir,jr]:=-h*a2*jcb[ir,jr];
j1[ir,ir]:=1-h*a2*jcb[ir,ir];
end;

{POLUCHAEM PRAVUYU CHAST'}
for ir:=1 to 5 do
y1[ir]:=y[ir]+b1*k1[ir];
dy(y1,f1);

for ir:=1 to 5 do
f1[ir]:=f1[ir]*h;
{RESHAEM SISTEMY}
kramer(j1,f1,k2);
{ --------- k2 NAIDENO ----------- }


{ ------- NAHODIM ynext ---------- }
for ir:=1 to 5 do
ynext[ir]:=y[ir]+w1*k1[ir]+w2*k2[ir];
xnext:=x+h;
end;


{-------------------}
begin
{ZAVODIM FAIL , V KOTORYI BUDET VYVODITSYA RESULTATY RASCHETOV}
assign(res,'c:\rsbtest1.dat');
rewrite(res);
{NACHALNYE USLOVIYA}
for i:=1 to 5 do
y[i]:=1;
x:=0;
h:=0.0005;
{ ----- RESHAEM SISTEMU ---- }
for i:=1 to 2000 do begin
rsnbrk(h,x,xnext,y,ynext,jcb);
{ZAPISYVAEM RESULTAT V FAIL,
ZAPISYVAEM TOLKO KAZHDUYU 10-YU TOCHKU}
if i mod 10 =0 then begin
write(res,x,'  ');
write(res,ynext[1],'  ');
write(res,ynext[2],'  ');
write(res,ynext[3],'  ');
write(res,ynext[4],'  ');
writeln(res,ynext[5]);
end;

for j:=1 to 5 do
y[j]:=ynext[j];
x:=xnext;
end;
{  ---- RESHENIE SISTEMY NAJDENO ----  }
close(res);
writeln;writeln;writeln;
writeln('                 SCHET ZAKONCHEN');
writeln;
write('          press < ENTER > to exit the program');readln;
end.

Немного теории можно посмотреть здесь:

http://www.library.cornell.edu/nr/bookcpdf/c16-6.pdf

http://www.ma.utexas.edu/documentation/nr/bookfpdf/f16-6.pdf

краткое описание использованного мной метода лежит здесь:
http://www.ualberta.ca/CMENG/courses/che674/handouts/chap6.pdf