Given the IVP

^{(1)}(

*t*) = f(

*t*, y(

*t*) )

y(

*a*) =

*y*

_{0}

where we want to estimated y(*b*) for *b* > *a*,
we can apply the same strategy we used for the composite-trapezoidal rule:

Break the interval [*a*, *b*] into *n* sub-intervals
and define *h* = (*b* - *a*)/*n*. Then set
*t*_{i} = *a* + *ih* for *i* = 0, 1, ..., *n*.
Therefore *t*_{0} = *a* and *t*_{n} = *b*.

Now, *y*_{0} is the initial value so let
*y*_{i} represent the approximation of
y(*t*_{i}) for *i* = 1, 2, ..., *n*.

Using any of the three techniques we've seen, Euler's, Heun's, or 4th-order Runge Kutta, we may now proceed as follows:

For *i* from 1 to *n*:

Approximate *y*_{i} using *y*_{i − 1}.

More specifically:

### Multiple-step Euler's Method

Calculate

*y*

_{i}=

*y*

_{i − 1}+

*h*f(

*t*

_{i − 1},

*y*

_{i − 1})

for *i* = 1, 2, ..., *n*.

### Multiple-step Heun's Method

Calculate

*K*_{0} = f(*t*_{i − 1}, *y*_{i − 1})

*K*_{1} = f(*t*_{i}, *y*_{i − 1} + *h* *K*_{0})

and set *y*_{i} = *y*_{i − 1} + *h* (*K*_{0} + *K*_{1})/2.

### 4th-order Runge Kutta

Calculate

*K*_{0} = f(*t*_{i − 1}, *y*_{i − 1})

*K*_{1} = f(*t*_{i − 1} + ½*h*, *y*_{i − 1} + ½*h* *K*_{0})

*K*_{2} = f(*t*_{i − 1} + ½*h*, *y*_{i − 1} + ½*h* *K*_{1})

*K*_{3} = f(*t*_{i}, *y*_{i − 1} + *h* *K*_{2})

and set *y*_{i} = *y*_{i − 1} + *h* (*K*_{0} + 2 *K*_{1} + 2 *K*_{2} + *K*_{3})/6.

Copyright ©2005 by Douglas Wilhelm Harder. All rights reserved.