Интеграл в правой части (8.4) вычисляется с помощью численной квадратуры

y_i^(k) – y₀ = h A_j⁽ⁱ⁾ f(x_j, y_j^(k-1)). (8.5)

Коэффициенты A_j⁽ⁱ⁾ системы (8.5) находятся из решения системы линейных алгебраических уравнений:

= A_j⁽ⁱ⁾ j^k-1; k=1, …n+1; i=1, …n. (8.6)

Система (8.6) получена из условия, что формула (8.5) точна для всех функций вида x, x², x³, …, xⁿ⁺¹ в точках x_i = ih, i=0, …, n.

Пусть n =2, дано три точки (x₀, x₁, x₂). Тогда матрица коэффициентов A_j⁽ⁱ⁾ равна

A= .

Для m =4, пять точек (x₀, x₁, x₂, x₃, x₄).

A= .

Численная реализация может быть выполнена следующим образом.

Шаг 1. Исходный интервал делим на 4 части, т. е. задаем m = 4. Определяем сетку x_i=a+ih₀, i =0, 1, …, 4, где h₀ = (a+b)/ 4 - шаг исходной сетки. Вычисляем по формулам (8.5) y_i⁽ⁿ⁾, i= 1, …, 4, пока не выполнится условие max < e, по всем i, где e - заданная точность.

Шаг 2. Делим интервал на 8 частей, т. е. уменьшаем шаг вдвое и находим y_i на сетках:

= a+ih₁, i =1, …, 4; h₁=h₀ /2;

= (b+a)/2+ih₁, i =1, …, 4.

Затем проверяем условие max < e в узлах первой сетки (с шагом h₀). Если оно выполняется, то вычисления заканчиваем. Если не выполняется, то исходный интервал уменьшаем вдвое, т. е. a₁=a; b₁=(a+b)/2, и переходим на шаг 1, где a=a₁; b=b₁.