Розглянемо звичайне диференціальне рівняння: dx/dt = f (t, x), у якому:
t — дійсна незалежна змінна (час);
x — змінна, залежна від t, що набуває величин у Rⁿ;
f — функція з R^n + 1 в Rⁿ, що задовольняє умови теореми існування і єдиності розв'язку задачі Коші.

При лінійній залежності f від x можна явно вирахувати розв'язки рівняння. У випадку нелінійної залежності при наближеному інтегруванні за допомогою різницевих схем бажано інтегрувати з якомога більшим кроком, щоб зменшити кількість обчислень нелінійного поля напрямків f. Відомі методи Рунге-Кутти різних порядків точності для інтегрування рівнянь. Похибка на кроці інтегрування h методу m-го порядку складає Ah^m + 1. При сталому кроці інтегрування без знання верхньої межі для сталої A неможливо вибрати максимальний крок інтегрування, при якому похибка інтегрування на одному кроці не перевищує певну сталу C.

Запровадимо такі позначення:

x_t — величина x певній величині t ;

y₁ = ¹/₃ h f (t, x_t) ;

y₂ = ¹/₃ h f (t + ¹/₃ h, x_t+ y₁) ;

y₃ = ¹/₃ h f (t + ¹/₃ h, x_t+ ¹/₂ y₁ + ¹/₂ y₂) ;

y₄ = ¹/₃ h f (t + ¹/₂ h, x_t + ³/₈ y₁ + ⁹/₈ y₃) ;

y₅ = ¹/₃ h f (t + h, x_t + ³/₂ y₁ – ⁹/₂ y₃ + 6y₄) ;

x_{t + h} = x_t + ¹/₂(y₁ + 4y₄ + y₅) + O(h⁵) — величина x для часу (t + h) ;

d = y₁ – ⁹/₂ y₃ + 4y₄ – ¹/₂ y₅ — похибка кроку інтегрування.

У численній літературі можна зустріти формули, що відрізняються сталими множниками для проміжних величин, але збігаються щодо кінцевого результату. Така модифікація Мерсона методу Рунге-Кутти (R.H. Merson, «An operational method for the study of integration processes», Proc. Symp. Data Processing , Weapons Res. Establ. Salisbury, Salisbury (1957) p. 110–125) передбачає автоматичний вибір кроку інтегрування:

• при справдженні нерівності || d || < ⁵/₃₂ C інтегрування продовжують з нової точки x_t + h з подвоєним кроком інтегрування 2h;

• при справдженні нерівностей ⁵/₃₂ C ≤ || d || ≤ 5C інтегрування продовжують з нової точки x_t + h з тим самим кроком інтегрування h;

• при справдженні нерівності 5C < || d || інтегрування продовжують з попередньої точки x_t зі зменшеним удвічі кроком інтегрування h/2.

Тут || d || — величина норми (довжини) вектора d = (d₁; d₂; d₃; …; d_n).

Відомі такі норми для векторних просторів R ⁿ :

|| d ||_p = ( | d₁| ^p + | d₂| ^p + | d₃| ^p + … + | d_n| ^p) ^1/p при 1 ≤ p.

Найменшу кількість обчислень здійснюють при обчисленні такої норми:

|| d ||₁ = | d₁| + | d₂| + | d₃| + … + | d_n|.

Подамо приклад програмної реалізації мовою Turbo Pascal 7.0.

{$N+}
const n=2;
type num=extended;
     arr=array[1..n] of num;
var
f,  {поле напрямків}
x,  {точка, у якiй вираховують f}
x0, {початкове значення, розв'язок}
y1,y2,y3,y4,y5: arr;    {див. опис}
t,               {час - аргумент f}
t0,           {поточний час для x0}
t1,           {кiнцевий час}
h,              {крок інтегрування}
err,   {похибка кроку інтегрування}
errup: num;       {верхня межа err}
d,g,nx,ny: integer;
j: byte;
o: text;            {вихiдний файл}

procedure flow; BEGIN
f[1]:=-x[2];
f[2]:= x[1];      END;

procedure merson;            BEGIN
             if (t0>t1) then begin
writeln('Wrong time limits!');
                          halt end;
repeat
 if t0+h>t1 then h:=t1-t0;
 for j:=1 to n do x[j]:=x0[j];
 flow;
 t:=t0+h/3;
 for j:=1 to n do            begin
 y1[j]:=f[j]*h/3;
  x[j]:=x0[j]+y1[j]            end;
 flow;
 for j:=1 to n do            begin
 y2[j]:=f[j]*h/3;
  x[j]:=x0[j]+(y1[j]+y2[j])/2  end;
 flow;
 t:=t0+h/2;
 for j:=1 to n do            begin
 y3[j]:=f[j]*h/3;
  x[j]:=x0[j]+3*y1[j]/8
             +9*y3[j]/8        end;
 flow;
 t:=t0+h;
 for j:=1 to n do            begin
 y4[j]:=f[j]*h/3;
  x[j]:=x0[j]+3*y1[j]/2-9*y3[j]/2
             +6*y4[j]          end;
 flow;
 err:=0;
 for J:=1 to n do            begin
 y5[j]:=f[j]*h/3;
 err:=err+abs(y1[j]-9*y3[j]/2
           +4*y4[j]-  y5[j]/2) end;
 if err>errup*5 then h:=h/2
                      else   begin
 for j:=1 to n do            begin
 x0[j]:=x0[j]+(y1[j]
      +4*y4[j]+y5[j])/2        end;
  t0:=t;
 if err < errup*0.15625
 then h:=h*2 end
until t0 > t1;                END;
BEGIN
errup:=1.E-13;
x0[1]:=1;
x0[2]:=0;
t0:=0;
t1:=33*pi;
h:=1.;
merson;
assign(o,'OUTPUT.txt');
append(o);
writeln(o,x0[1],' ',x0[2],' ',
        sqr(x0[1])+sqr(x0[2]));
close(o);
END.

За поле напрямів вибрано поле поворотів навколо початку координат, що зберігає суму квадратів динамічних змінних.

Для типу змінних real маємо такий рядок вихідного файлу.

-1.00000000000364E+0000  5.70762659357626E-0010  1.00000000000728E+0000

Для типу змінних extended маємо наступний рядок вихідного файлу.

-1.00000000000000E+0000  5.36411451727628E-0010  1.00000000000000E+0000