Chapter 8 · Section 8.8

8.8 Numerical Methods V: Higher-Order Equations and Systems

The methods of Sections 8.4–8.7 were written for a single first-order equation $y^{'} = f (x, y)$ . Engineering models are rarely that simple. A mass on a spring, an RLC circuit or a swinging pendulum obeys a second-order equation, and a structure with several masses or a circuit with several loops leads to a system of equations. Most of these problems have no formula solution, especially when they are nonlinear. The good news is that nothing new has to be invented: every higher-order equation can be rewritten as a first-order system, and every method we know works for systems with essentially no change. This section shows how.

Because the applications are usually about motion in time, from now on we use $t$ for the independent variable and write the unknowns as $x (t)$ , $y (t)$ , as in Chapter 7.

From one higher-order equation to a first-order system

Consider a second-order equation solved for the highest derivative,

x^{''} = F (t, x, x^{'}), x (t_{0}) = α, x^{'} (t_{0}) = β .

Introduce the velocity as a new unknown, $y = x^{'}$ . Then $y^{'} = x^{''}$ , and the single equation becomes the pair

x^{'} y^{'} = y, = F (t, x, y), x (t_{0}) y (t_{0}) = α, = β .

The two problems are equivalent: if $x$ solves the second-order problem, then $(x, x^{'})$ solves the system; conversely, the first equation of the system says $y = x^{'}$ , and then the second says $x^{''} = F (t, x, x^{'})$ . The two initial conditions of the original problem turn into one initial condition for each unknown.

Definition 8.8.1 (Equivalent first-order system)

For the $n$ th-order equation $x^{(n)} = F (t, x, x^{'}, \dots, x^{(n - 1)})$ with $x^{(k)} (t_{0}) = α_{k}$ ( $k = 0, \dots, n - 1$ ), set

u_{1} = x, u_{2} = x^{'}, \dots, u_{n} = x^{(n - 1)} .

Then

u_{1}^{'} u_{2}^{'} u_{n - 1}^{'} u_{n}^{'} = u_{2}, = u_{3}, ⋮ = u_{n}, = F (t, u_{1}, u_{2}, \dots, u_{n}), u_{k} (t_{0}) = α_{k - 1} .

In vector form, $u^{'} = F (t, u)$ , $u (t_{0}) = u_{0}$ .

Example 8.8.1 (A third-order equation)

Write the initial-value problem $x^{'''} + t x^{''} - 2 (x^{'})^{2} + x = sin t, x (0) = 1, x^{'} (0) = 0, x^{''} (0) = - 2,$ as a first-order system.

Solution

Solve for the highest derivative. $x^{'''} = sin t - t x^{''} + 2 (x^{'})^{2} - x$ .
New unknowns. $u_{1} = x$ , $u_{2} = x^{'}$ , $u_{3} = x^{''}$ .
The system. $u_{1}^{'} u_{2}^{'} u_{3}^{'} = u_{2}, = u_{3}, = sin t - t u_{3} + 2 u_{2}^{2} - u_{1}, u_{1} (0) u_{2} (0) u_{3} (0) = 1, = 0, = - 2.$ Note that the system is nonlinear (because of $u_{2}^{2}$ ); that is no obstacle for numerical methods.

Systems of several higher-order equations are handled the same way: introduce one new unknown for each derivative below the highest order of each unknown function (Exercise 2).

Numerical methods for systems

We write the general two-equation system as

x^{'} = f (t, x, y), y^{'} = g (t, x, y), x (t_{0}) = x_{0}, y (t_{0}) = y_{0} .

The idea is simply to treat the pair $(x, y)$ as one unknown vector. Every formula for $y^{'} = f (x, y)$ is applied componentwise, with the important rule that all components are advanced together, using values from the same stage.

Euler's method

x_{n + 1} y_{n + 1} = x_{n} + h f (t_{n}, x_{n}, y_{n}), = y_{n} + h g (t_{n}, x_{n}, y_{n}) .

Both right-hand sides use the old values $(x_{n}, y_{n})$ .

The RK4 method

For each stage we compute a slope $k_{i}$ for $x$ and a slope $l_{i}$ for $y$ , and both are evaluated at the same trial point:

k_{1} k_{2} k_{3} k_{4} = f (t_{n}, x_{n}, y_{n}), = f (t_{n} + \frac{h}{2}, x_{n} + \frac{h}{2} k_{1}, y_{n} + \frac{h}{2} l_{1}), = f (t_{n} + \frac{h}{2}, x_{n} + \frac{h}{2} k_{2}, y_{n} + \frac{h}{2} l_{2}), = f (t_{n} + h, x_{n} + h k_{3}, y_{n} + h l_{3}), l_{1} l_{2} l_{3} l_{4} = g (t_{n}, x_{n}, y_{n}), = g (t_{n} + \frac{h}{2}, x_{n} + \frac{h}{2} k_{1}, y_{n} + \frac{h}{2} l_{1}), = g (t_{n} + \frac{h}{2}, x_{n} + \frac{h}{2} k_{2}, y_{n} + \frac{h}{2} l_{2}), = g (t_{n} + h, x_{n} + h k_{3}, y_{n} + h l_{3}),

and then

x_{n + 1} y_{n + 1} = x_{n} + \frac{h}{6} (k_{1} + 2 k_{2} + 2 k_{3} + k_{4}), = y_{n} + \frac{h}{6} (l_{1} + 2 l_{2} + 2 l_{3} + l_{4}) .

In vector notation this is exactly Definition 8.6.1 with $u = (x, y)$ , $F = (f, g)$ and vector slopes $k_{i} = (k_{i}, l_{i})$ . The same is true of Heun's method and of ABAM, and the orders do not change: Euler is still first order, Heun second, RK4 and ABAM fourth.

Method in a nutshell: solving a higher-order IVP numerically

Solve the equation for its highest derivative.
Introduce $x$ , $y = x^{'}$ (and $z = x^{''}$ , … for higher order); write the first-order system and its initial values.
Choose a method and a step $h$ ; check that $h$ is small compared with the time scales of the problem (see the section on stability below).
At every stage compute all component slopes at the same trial point, then advance all components together.
Read off $x_{n}$ (the original unknown) and, as a bonus, $y_{n} \approx x^{'} (t_{n})$ .

Common mistake

Do not update $x$ first and then use the new $x_{n + 1}$ while computing $y_{n + 1}$ (or the new $k_{2}$ when computing $l_{2}$ ). Within a step, every slope of a given stage must be computed from the same point. Mixing old and new values silently produces a different method, usually with a lower order.

Worked example: a damped oscillator

A mass of $1$ kg hangs on a spring of stiffness $10$ N/m with a damper of coefficient $2$ N·s/m. It is pulled down $1$ m and released from rest. Its displacement satisfies (Chapter 5)

x^{''} + 2 x^{'} + 10 x = 0, x (0) = 1, x^{'} (0) = 0.

This linear problem has a known solution, which lets us measure the errors. The characteristic equation $r^{2} + 2 r + 10 = 0$ has roots $r = - 1 \pm 3 i$ , so $x = e^{- t} (c_{1} cos 3 t + c_{2} sin 3 t)$ . The initial conditions give $c_{1} = 1$ and $- c_{1} + 3 c_{2} = 0$ , so

x (t) = e^{- t} (cos 3 t + \frac{1}{3} sin 3 t), x^{'} (t) = - \frac{10}{3} e^{- t} sin 3 t .

With $y = x^{'}$ the equivalent system is

x^{'} = y, y^{'} = - 10 x - 2 y, x (0) = 1, y (0) = 0.

So $f (t, x, y) = y$ and $g (t, x, y) = - 10 x - 2 y$ .

Example 8.8.2 (Euler's method for the oscillator)

Take two Euler steps with $h = 0.1$ and compare with the exact values.

Solution

First step from $(x_{0}, y_{0}) = (1, 0)$ : $x_{1} y_{1} = 1 + 0.1 (0) = 1, = 0 + 0.1 (- 10 \cdot 1 - 2 \cdot 0) = - 1.$
Second step from $(x_{1}, y_{1}) = (1, - 1)$ : $x_{2} y_{2} = 1 + 0.1 (- 1) = 0.9, = - 1 + 0.1 (- 10 \cdot 1 - 2 (- 1)) = - 1 + 0.1 (- 8) = - 1.8.$
Compare. Exact: $x (0.1) = 0.953557$ , $x^{'} (0.1) = - 0.891326$ ; $x (0.2) = 0.829824$ , $x^{'} (0.2) = - 1.540967$ . After two steps the position is already off by $0.07$ . Notice that the first step leaves $x$ unchanged because the velocity starts at zero: Euler cannot "see" that the mass begins to accelerate.

Example 8.8.3 (RK4 for the oscillator)

Solve the same problem by RK4 with $h = 0.1$ on $0 \leq t \leq 1$ . Show the first step in full.

Solution

Stage 1 at $(t, x, y) = (0, 1, 0)$ : $k_{1} = y = 0, l_{1} = - 10 (1) - 2 (0) = - 10.$
Stage 2 at $(0.05, 1 + 0.05 (0), 0 + 0.05 (- 10)) = (0.05, 1, - 0.5)$ : $k_{2} = - 0.5, l_{2} = - 10 (1) - 2 (- 0.5) = - 9.$
Stage 3 at $(0.05, 1 + 0.05 (- 0.5), 0 + 0.05 (- 9)) = (0.05, 0.975, - 0.45)$ : $k_{3} = - 0.45, l_{3} = - 10 (0.975) - 2 (- 0.45) = - 8.85.$
Stage 4 at $(0.1, 1 + 0.1 (- 0.45), 0 + 0.1 (- 8.85)) = (0.1, 0.955, - 0.885)$ : $k_{4} = - 0.885, l_{4} = - 10 (0.955) - 2 (- 0.885) = - 7.78.$
Combine. $x_{1} y_{1} = 1 + \frac{0.1}{6} (0 + 2 (- 0.5) + 2 (- 0.45) - 0.885) = 1 + \frac{0.1}{6} (- 2.785) = 0.953583, = 0 + \frac{0.1}{6} (- 10 + 2 (- 9) + 2 (- 8.85) - 7.78) = \frac{0.1}{6} (- 53.48) = - 0.891333.$ Compare with the exact $0.953557$ and $- 0.891326$ : errors of $2.7 \times 1 0^{- 5}$ and $7.5 \times 1 0^{- 6}$ , against $0.046$ and $0.11$ for Euler.

All ten steps.

$t_{n}$	$x_{n}$	exact $x (t_{n})$	$y_{n}$	exact $x^{'} (t_{n})$
0.1	0.953583	0.953557	−0.891333	−0.891326
0.2	0.829874	0.829824	−1.541026	−1.540967
0.3	0.653997	0.653934	−1.934477	−1.934343
0.4	0.451214	0.451150	−2.082762	−2.082548
0.5	0.244627	0.244575	−2.016982	−2.016704
0.6	0.053492	0.053462	−1.781844	−1.781530
0.7	−0.107813	−0.107813	−1.429172	−1.428857
0.8	−0.230195	−0.230164	−1.011964	−1.011684
0.9	−0.309710	−0.309648	−0.579411	−0.579199
1.0	−0.346979	−0.346893	−0.173172	−0.173050

Over the whole interval the position error stays below

9 \times 1 0^{- 5}

and the velocity error below

4 \times 1 0^{- 4}

Figure 8.8.1. RK4 with

h = 0.1

for

x^{''} + 2 x^{'} + 10 x = 0

x (0) = 1

x^{'} (0) = 0

. Left: position and velocity against time. Right: the phase plane

(x, x^{'})

; the trajectory spirals into the equilibrium at the origin. The RK4 points lie on the exact curves to plotting accuracy.

Use the solver to reproduce Example 8.8.3 and to explore. Switch the method to Euler and try $h = 0.1$ and then $h = 0.25$ ; or change $g$ to model a different spring, for instance -25*x - 0.5*y.

A nonlinear example: the pendulum

For a linear oscillator we did not really need a computer. The power of numerical methods shows when the equation is nonlinear. A pendulum of length $L$ swinging through angle $θ (t)$ in air obeys

θ^{''} + c θ^{'} + \frac{g}{L} sin θ = 0.

For small swings one replaces $sin θ$ by $θ$ and obtains a linear equation like the one above. For large swings this approximation is poor, and the nonlinear equation has no elementary solution.

Example 8.8.4 (Large swings of a pendulum)

Take $g / L = 9 s^{- 2}$ (a pendulum about $1.09$ m long) and $c = 0.3 s^{- 1}$ . The pendulum is released from rest at $θ (0) = 1$ rad (about $5 7^{\circ}$ ). Use RK4 with $h = 0.05$ to compute $θ (t)$ for $0 \leq t \leq 3$ , and compare with the linearized model $θ^{''} + 0.3 θ^{'} + 9 θ = 0$ .

Solution

System. With $x = θ$ , $y = θ^{'}$ : $x^{'} = y, y^{'} = - 9 sin x - 0.3 y, x (0) = 1, y (0) = 0.$
First stage, as a check. $k_{1} = 0$ , $l_{1} = - 9 sin 1 = - 7.573239$ . The remaining stages follow the pattern of Example 8.8.3; the first step gives $x_{1} = 0.990590$ , $y_{1} = - 0.375070$ .
Linearized model. Its characteristic roots are $- 0.15 \pm ωi$ with $ω = 9 - 0.0225 = 2.996248$ , and with the same initial data $θ_{lin} (t) = e^{- 0.15 t} (cos ω t + \frac{0.15}{ω} sin ω t) .$

Results (60 RK4 steps; values every

0.5

s, rounded to four decimals):

$t$	$θ$ (RK4)	$θ^{'}$ (RK4)	$θ_{lin}$
0.5	0.2055	−2.6357	0.1137
1.0	−0.8010	−0.7271	−0.8454
1.5	−0.3850	2.0657	−0.2117
2.0	0.5951	1.1424	0.6991
2.5	0.4639	−1.5130	0.2764
3.0	−0.4117	−1.3263	−0.5645

A computation with a much smaller step agrees with these RK4 values to within

5 \times 1 0^{- 5}

Interpretation. The two models differ by more than $0.2$ rad (about $0.23$ rad at worst) already within the first three seconds. The nonlinear pendulum swings more slowly: for a swing of $1$ rad its period is about $7%$ longer than the small-angle value $2 π /3$ s, so it falls further and further behind the linear model (Figure 8.8.2).

Figure 8.8.2. Example 8.8.4. The nonlinear pendulum (solid, RK4) has a longer period than the linearized model (dashed), so the two curves drift out of phase.

Try the pendulum in the solver above: enter $g$ as -9*sin(x) - 0.3*y, $x_{0} = 1$ , $h = 0.05$ and $n = 60$ . Then raise $x_{0}$ to $3$ (almost upside down) and see how the motion changes.

Stability: when the step is too large

Example 8.4.6 gave a first warning: a larger $h$ does not always mean just a larger error. For oscillating or decaying solutions something more dramatic can happen: beyond a certain step size the numerical solution grows without bound, even though the true solution dies out. This is called numerical instability, and for systems it needs a closer look.

The test equation

The cleanest way to understand it is the test equation $y^{'} = λ y$ with $Re λ < 0$ , whose exact solution $y_{0} e^{λ t}$ decays. One step of a method multiplies $y_{n}$ by a fixed factor, which we found in Sections 8.4 to 8.6. With $z = hλ$ :

Euler: RK4: y_{n + 1} = (1 + z) y_{n}, y_{n + 1} = (1 + z + \frac{z ^{2}}{2} + \frac{z ^{3}}{6} + \frac{z ^{4}}{24}) y_{n} .

The numbers decay only if the factor has absolute value at most $1$ . For real $λ < 0$ this means

Euler: h < \frac{2}{∣ λ ∣}, RK4: h \leq \frac{2.785}{∣ λ ∣} (approximately) .

Why does this matter for systems? For a linear system $u^{'} = A u$ the solution is built from modes $e^{λ t}$ , one for each eigenvalue $λ$ of $A$ (Chapter 7), and a numerical method treats every mode like a separate test equation. So each eigenvalue must satisfy the stability condition. For a nonlinear system the same reasoning applies, approximately, to the eigenvalues of the linearization near the solution.

Example 8.8.5 (Euler's method blows up)

For the oscillator $x^{'} = y$ , $y^{'} = - 10 x - 2 y$ of Example 8.8.2, find the largest step for which Euler's method gives decaying numbers. What happens with $h = 0.25$ ?

Solution

Eigenvalues. The matrix $A = (0 - 10 1 - 2)$ has characteristic equation $λ^{2} + 2 λ + 10 = 0$ , so $λ = - 1 \pm 3 i$ (the same numbers as the characteristic roots of the original equation).
Euler's factor. For $λ = - 1 + 3 i$ : $∣1 + hλ ∣^{2} = ∣ (1 - h) + 3 h i ∣^{2} = (1 - h)^{2} + 9 h^{2} = 1 - 2 h + 10 h^{2} .$ (The conjugate eigenvalue gives the same value.)
Condition. $1 - 2 h + 10 h^{2} < 1 ⟺ 10 h^{2} - 2 h < 0 ⟺ 0 < h < 0.2$ .
With $h = 0.25$ . The factor is $∣0.75 + 0.75 i ∣ = 1.0607$ per step. After $32$ steps ( $t = 8$ ) Euler gives $x_{32} = 6.58$ , while the true displacement is only about $4 \times 1 0^{- 5}$ . With $h = 0.1$ the factor is $∣0.9 + 0.3 i ∣ = 0.949$ and the numbers decay, although more slowly than the true factor $e^{- 0.1} = 0.905$ .
RK4. Solving $∣ R (hλ) ∣ = 1$ numerically for $λ = - 1 + 3 i$ gives $h \approx 0.89$ . RK4 stays stable with steps more than four times larger than Euler can use.

Figure 8.8.3. Left: Euler's method for

x^{''} + 2 x^{'} + 10 x = 0

decays with

h = 0.1

but grows with

h = 0.25 > 0.2

. Right: for the undamped oscillator

x^{''} + 4 x = 0

, Euler's points spiral outward for every

h

, while RK4 stays on the exact closed orbit.

Undamped oscillations

For $x^{''} + ω^{2} x = 0$ the eigenvalues are $\pm ωi$ and Euler's factor is $∣1 + iωh ∣ = 1 + ω^{2} h^{2} > 1$ for every $h > 0$ : Euler always pumps energy into an undamped oscillator (Exercise 5). RK4's factor satisfies $∣ R (iωh) ∣^{2} = 1 - \frac{( ωh ) ^{6}}{72} + \frac{( ωh ) ^{8}}{576}$ , which is at most $1$ when $ωh \leq 22$ ; RK4 loses a tiny amount of energy per step instead. This is one reason why Euler's method is never used for long simulations of vibrating systems.

Stiff problems

Trouble arises when a system contains time scales of very different sizes, for example eigenvalues $- 1$ and $- 1000$ . The fast mode dies out almost immediately and is invisible in the solution, yet an explicit method must still respect $h ≲ 2.8/1000$ for RK4, for the whole computation. Such problems are called stiff. They are common in chemical kinetics and circuits, and they are solved with implicit methods (relatives of the implicit trapezoid and Adams–Moulton formulas), which are beyond this course. A practical warning sign: if halving $h$ changes the answer completely, or the numbers explode, suspect instability before blaming the model.

Summary

An $n$ th-order equation becomes a system of $n$ first-order equations by taking $x, x^{'}, \dots, x^{(n - 1)}$ as unknowns.
Euler, Heun, RK4 and ABAM work for systems componentwise; all components are advanced together, stage by stage. Their orders are unchanged.
For decaying or oscillating problems the step must also satisfy a stability condition, $∣ R (hλ) ∣ \leq 1$ for every eigenvalue $λ$ . Euler requires $h < 2/∣ λ ∣$ for real $λ < 0$ ; RK4 allows about $2.785/∣ λ ∣$ .
Euler's method is unstable for undamped oscillations at every step size.

Exercises

Write $x^{'''} - t x^{'} + x^{2} = cos t$ , $x (0) = 1$ , $x^{'} (0) = 0$ , $x^{''} (0) = 2$ as a first-order system.
Answer
$u_{1}^{'} = u_{2}$ , $u_{2}^{'} = u_{3}$ , $u_{3}^{'} = cos t + t u_{2} - u_{1}^{2}$ ; $u_{1} (0) = 1$ , $u_{2} (0) = 0$ , $u_{3} (0) = 2$ .
Two coupled masses obey $x_{1}^{''} = - 5 x_{1} + 2 x_{2}$ , $x_{2}^{''} = 2 x_{1} - 2 x_{2}$ . Write this as a system of four first-order equations.
Answer
With $u_{1} = x_{1}$ , $u_{2} = x_{1}^{'}$ , $u_{3} = x_{2}$ , $u_{4} = x_{2}^{'}$ : $u_{1}^{'} = u_{2}$ , $u_{2}^{'} = - 5 u_{1} + 2 u_{3}$ , $u_{3}^{'} = u_{4}$ , $u_{4}^{'} = 2 u_{1} - 2 u_{3}$ .
Take two Euler steps with $h = 0.1$ for $x^{''} + 4 x = 0$ , $x (0) = 1$ , $x^{'} (0) = 0$ . Compare with $x = cos 2 t$ .
Answer
$(x_{1}, y_{1}) = (1, - 0.4)$ , $(x_{2}, y_{2}) = (0.96, - 0.8)$ ; exact $x (0.2) = 0.921061$ , $x^{'} (0.2) = - 0.778837$ .
Take one RK4 step with $h = 0.5$ for $x^{'} = y$ , $y^{'} = - x$ , $x (0) = 1$ , $y (0) = 0$ . Show that $x_{1} = 1 - \frac{h ^{2}}{2} + \frac{h ^{4}}{24}$ and $y_{1} = - h + \frac{h ^{3}}{6}$ , and compare with $cos 0.5$ and $- sin 0.5$ .
Answer
$x_{1} = 0.877604$ ( $cos 0.5 = 0.877583$ ), $y_{1} = - 0.479167$ ( $- sin 0.5 = - 0.479426$ ). These are the Taylor polynomials of $cos h$ and $- sin h$ .
For $x^{'} = y$ , $y^{'} = - ω^{2} x$ , show that Euler's method satisfies $ω^{2} x_{n + 1}^{2} + y_{n + 1}^{2} = (1 + ω^{2} h^{2}) (ω^{2} x_{n}^{2} + y_{n}^{2})$ . What does this mean physically?
Answer
$x_{n + 1} = x_{n} + h y_{n}$ , $y_{n + 1} = y_{n} - h ω^{2} x_{n}$ ; expanding, the cross terms cancel and the identity follows. The (scaled) energy grows by the factor $1 + ω^{2} h^{2}$ every step, so the computed oscillation spirals outward.
For $x^{''} + 4 x^{'} + 13 x = 0$ , find the largest step for which Euler's method gives decaying numbers.
Answer
$λ = - 2 \pm 3 i$ ; $(1 - 2 h)^{2} + 9 h^{2} < 1 ⟺ 13 h^{2} < 4 h ⟺ h < 4/13 \approx 0.308$ .
For $y^{'} = - 50 y$ , what is the largest stable step for Euler's method and (approximately) for RK4?
Answer
Euler: $h < 2/50 = 0.04$ . RK4: $h \leq 2.785/50 \approx 0.0557$ .
A predator–prey model reads $x^{'} = x - x y$ , $y^{'} = - y + x y$ with $x (0) = 2$ , $y (0) = 1$ . Take one RK4 step with $h = 0.1$ .
Answer
$(k_{1}, l_{1}) = (0, 1)$ , $(k_{2}, l_{2}) = (- 0.1, 1.05)$ , $(k_{3}, l_{3}) = (- 0.104738, 1.047238)$ , $(k_{4}, l_{4}) = (- 0.208351, 1.093153)$ ; $x_{1} = 1.989703$ , $y_{1} = 1.104794$ .
The charge in a series RLC circuit satisfies $q^{''} + 6 q^{'} + 25 q = 50$ , $q (0) = 0$ , $q^{'} (0) = 0$ . Write it as a system for $q$ and the current $i = q^{'}$ , and take one RK4 step with $h = 0.1$ . Compare with the exact solution $q = 2 - e^{- 3 t} (2 cos 4 t + 1.5 sin 4 t)$ , $i = 12.5 e^{- 3 t} sin 4 t$ .
Answer
$q^{'} = i$ , $i^{'} = 50 - 25 q - 6 i$ . $q_{1} = 0.202292$ , $i_{1} = 3.609167$ ; exact $q (0.1) = 0.202590$ , $i (0.1) = 3.606103$ .
Write Heun's method for the system $x^{'} = f (t, x, y)$ , $y^{'} = g (t, x, y)$ , and take one step with $h = 0.1$ for the oscillator of Example 8.8.2.
Answer
Predictor $(x^{*}, y^{*}) = (x_{n} + h k_{1}, y_{n} + h l_{1})$ , slopes $(k_{2}, l_{2})$ at $(t_{n} + h, x^{*}, y^{*})$ , then $x_{n + 1} = x_{n} + \frac{h}{2} (k_{1} + k_{2})$ , $y_{n + 1} = y_{n} + \frac{h}{2} (l_{1} + l_{2})$ . Here $(k_{1}, l_{1}) = (0, - 10)$ , $(x^{*}, y^{*}) = (1, - 1)$ , $(k_{2}, l_{2}) = (- 1, - 8)$ , so $(x_{1}, y_{1}) = (0.95, - 0.9)$ (exact $0.953557$ , $- 0.891326$ ).