# Preface

This section gives an introduction to polynomial approximations based on Taylor polynomials. In principle, it is possible to not only determine the numerical approximations for the values of the exact solution at the mesh points x0, x1, x2, … , but also find polynomial approximations in the intervals between discrete mesh points. By joining grid points (x0, y0), (x1, y1), … with straight-line segments, we are able to obtain an overall approximation to the solution curve corresponding to the solution of the given initial-value problem.

The Taylor series method is of general applicability, and it is the standard to which the accuracy of various other numerical methods are compared. It can be devised to have any specified degree of accuracy under some general assumptions. The most important single result in numerical computations is Taylor's theorem, which we now state below.

Theorem (Taylor): : Let f(x) have n+1 continuous derivatives on the closed interval [𝑎.b] for some positive integer n, and let $$x, x_0 \in [a,b] .$$ Then

$f(x) = p_n + R_n (x)$
for
$p_n (x) = f(x_0 ) + f' (x_0 ) \,(x-x_0 ) + \cdots + \frac{f^{(n)} (x_0 )}{n!} \, (x- x_0 )^n = \sum_{k=0}^n \frac{(x- x_0 )^k}{k!} \, f^{(k)} (x_0 )$
and
$R_n (x) = \frac{1}{n!} \int_{x_0}^x (x-t)^n (t)\, {\text d}t = \frac{(x- x_0 )^{n+1}}{(n+1)!} \,f^{(n+1)} (\xi_x ) ,$
where ξ is a point between x and x0.     ⧫

Brook Taylor (1685--1731) was educated at St. John's College of the Cambridge University (England), from which he graduated in 1709. He wrote two very significant books, ‘Methodus incrementorum directa et inversa’ and ‘Linear Perspective,’ which were published in 1715. He was the individual to invent ‘Integration of Parts’ and also a series called the Taylor’s Expansion. The Taylor’s Theorem is based on the letter written by him to Machin in which he told about the origination of this idea. However, it appears that he did not appreciate its importance and he certainly did not bother with a proof of Taylor's theorem, which was discovered and published in 1694 by Johann Bernoulli (1667--1748). Taylor acknowledged that his work was based on that of Newton and Kepler. However, Johann Bernoulli claimed his priority in discovering both of his main results (Taylor's expansion and integration by parts).

 Johann Bernoulli Brook Taylor

# Polynomial Approximations

The Taylor series method is of general applicability, and it is a standard to which we compare the accuracy of the various other numerical methods for solving an initial value problem for ordinary differential equations:
$y' = f(x,y ) , \quad y(x_0 ) = y_0 .$
A numerical method for solving the above initial value problem starts with choosing a set of mesh points
$\left( x_0 , y_0 \right) , \left( x_1 , y_1 \right) , \ldots , \left( x_N , y_N \right) ,$
where $$x_n = x_0 + n\,h ,$$ n = 0, 1, …. For simplicity, we consider the set of uniformly distributed grid points starting from x0. We shall always denote by yn a number intended to approximate y(xn), the value of the exact solution of the initial value problem y' = f(x,y), y(x0) = y0 at x = xn. In what follows, it will be assumed that the given slope function f(x,y) has as many derivatives in a neighborhood of the initial point as needed.

Theorem: Suppose that a function y(x) is a solution of the initial value problem $$y' = f(x,y), \ y(x_0 ) = y_0$$ and that it has m+1 continuous derivatives on the interval [a,b] containing grid points $$\{ x_n \}_{n\ge 0} .$$ Assume that y(x) has a Taylor series expansion of order m about fixed value xn ∈ [𝑎,b] :

$y(x_n +h) = y(x_n ) + h\, \Phi_m \left( x_n, y(x_n )\right) + O\left( h^{m+1} \right) ,$
where the increment function is
$\Phi_m \left( x_n, y(x_n )\right) = \sum_{j=1}^m \frac{y^{(j)} (x_n )}{j!} \, h^{j-1}$
and
$y^{(k)} (x) = \left( \frac{\partial}{\partial x} + f\,\frac{\partial}{\partial y} \right)^{k-1}\, f\left( x, y(x) \right) , \qquad k=0,1,2,\ldots .$
In particular,
\begin{align*} y' (x) &= f(x,y(x)) , \\ y'' (x) &= f_x + f_y \, y' = f_x + f\, f_y , \\ y^{(3)} (x)&= f_{xx} + 2\, f_{xy} y' + f_y y'' + f_{yy} \left( y' \right)^2 \\ &= f_{xx} + 2\,f_{xy}\,f + f_{yy} f^2 + f_y \left( f_x + f\, f_y \right) , \\ y^{(4)} (x)&= f_{xxx} + 3\, f_{xxy} y' +3\, f_{xyy} \left( y' \right)^2 +3\, f_{xy} y'' + f_y y''' +3\, f_{yy} y' y'' + f_{yyy} \left( y' \right)^3 \\ &= f_{xxx} + 3\, f_{xxy} \,f +3\, f_{xyy} \left( f \right)^2 +3\, f_{xy} \left( f_x + f\, f_y \right) + f_{yy} f\left( f_x + f\, f_y \right) + f_y y''' . \end{align*}
Then the approximate numerical solution at mesh points $$\left( x_n , y_n \right) , \ n=0,1,2,\ldots , N,$$ to the initial value problem $$y' = f(x,y), \ y(x_0 ) = y_0$$ can be evaluated according to the formula:
$y_{n+1} = y_n + h \,f' (x_n ) + \frac{h^2}{2!}\, f'' (x_n ) + \cdots + \frac{h^m}{m!} \,f^{(m)} (x_n ) ,$
where h is a step size. This scheme is known as the Taylor algorithm of order m.     ⧫

The approximate numerical solution to the initial value problem $$y' = f(x,y ) , \quad y(x_0 ) = y_0$$ over the interval [𝑎,b] is derived by using the Taylor series expansion on each subinterval $$[ x_n , x_{n+1}] .$$ The general step for Taylor's method of order m is

$y_{n+1} = y_n + h\,y' (x_n ) + \frac{h^2}{2} \,y'' (x_n ) + \cdots + \frac{h^m}{m!} \, y^{(m)} (x_n ).$
The Taylor method of order m has the property that the final global error is of the order $$O\left( h^{m+1} \right) ,$$ that is,
$\left\vert y_n y (x_n ) \right\vert \le \frac{h^m}{(m+1)!}\, x_n e^{x_n} .$
This relation shows that at a fixed point xn, the error in the approximation of the exact solution y by the Taylor algorithm of order m tends to zero like hm as the stepsize h tends to zero.
It is possible theoretically to choose m as large as necessary to make this error as small as desired. However, in practice, we usually compute two sets of approximations using step sizes h and h/2 and compare the results.

Theorem (Precision of Taylor's method of order m): Suppose that the solution to the initial value problem $$y' = f(x,y ) , \quad y(x_0 ) = y_0$$ possesses m+1 continuous derivatives. Let $$\left( x_n , f_n \right) , \ n=0,1,2,\ldots , N,$$ be the sequence of approximations generated by Taylor's method of order m. Then

$\begin{split} \left\vert e_n \right\vert &= \left\vert y(x_n ) - y_n \right\vert = O \left( h^{m+1} \right) , \\ \left\vert e_{n+1} \right\vert &= \left\vert y(x_{n+1} ) - y_n - h\, \Phi_m \left( x_n, y(x_n )\right) \right\vert = O \left( h^{m} \right) . \end{split}$
In particular, the final global error at the end interval will satisfy
$\left\vert y(b) - e_m \right\vert = O \left( h^{m} \right) . \qquad \blacksquare$

Second Order Polynomial Approximations

Since the first order Taylor series approximation is identical with Euler’s method, we start with the second order one:

$y_{n+1} = y_n + h\,f (x_n , y_n ) + \frac{h^2}{2} \left[ f_x (x_n , y_n ) + f(x_n , y_n )\, f_y (x_n , y_n ) \right] = y_n + h\, \Phi_2 (h) ,$
where the increment function $$\displaystyle \Phi_2 (h) = \frac{h}{2} \left. \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right) f(x,y) \right\vert_{x=x_n , y= y_n}$$ is just adding the second order differential deviation to the next term in the same equation used for the first order Taylor expansion. We show how it works in the following examples.
Example: We start with the initial value problem for the Riccati equation
$y' = y^2 - x^2 , \qquad y(0) = 1/2.$
Upon taking the mesh of points with uniform step size h = 0.1, $$x_n = n\,h , \quad n=0,1,2,\ldots ,$$ we denote by yn approximate values of the unknown solution at these grid points. The initial point is (0, 0.5) is obtained from the given initial condition. For a second order Taylor's approximation, we only need to determine the first derivative of the slope function:
\begin{align*} f' (x,y) &= \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right) f(x,y) = 2y \left( y^2 - x^2 \right) - 2x = 2\,y^3 - 2\,y\,x^2 - 2\,x , \end{align*}
With this formula, we are able to make approximating increments: $y_{n+1} = y_n + h \left( y^2_n - x^2_n \right) + h^2 \left( y_n^3 - x_n^2 y_n - x_n \right) .$ Now we ask Mathematica to perform calculations:
Clear[y]
y[0] = 1/2; h = 0.1; f[x_, y_] = y^2 - x^2;
Do[y[n + 1] = y[n] + h f[h n, y[n]] + h^2 *(y[n]^3 - n*h - n^2 *h^2 *y[n]), {n, 0, 10}]
y[5]
which gives the Taylor polynomial approximation at point x = 0.5 to be y5 = 0.619956, which is closed to the true value of 0.6177691815444183.
The results of applying Taylor's approximation of order 2 with h = 0.1
n xn yn Exact solution Absolute error
1 0.1 0.52625 0.525974 -0.000276406
2 0.2 0.553349 0.552738 -0.000610995
3 0.3 0.579441 0.578417 -0.00102377
4 0.4 0.60244 0.600901 -0.0015396
5 0.5 0.619956 0.617769 -0.00218714
6 0.6 0.629224 0.626228 -0.00299599
7 0.7 0.627042 0.623052 -0.00399052
8 0.8 0.609753 0.604575 -0.00517862
9 0.9 0.573298 0.566763 -0.00653424
10 1.0 0.513405 0.505431 -0.00797472
■
Example: Consider the initial value problem $$y’=1/(2x-3y+5), \quad y(0)=1$$ and use the standard mathematica command DSolve. However, mathematica fails to provide the explicit solution because it is expressed through a special function, called the Lambert function (which we denote by W and which is defined from the equation $$W\,e^W = x$$ ):
$y(x) = \frac{2x}{3} + \frac{7}{8} - \frac{1}{2} \, W \left( \frac{1}{3}\,e^{(1+4x)/3} \right) .$
Let us try to solve this problem numerically using Taylor's algorithm of order 2. However, we check with a standard (very reliable) mathematica procedure, called NDSolve that provide a numerical solution to the given problem:
sol = NDSolve[{Derivative[1][y][x] == 1/(2*x - 3*y[x] + 5), y[0] == 1}, y, {x, 0, 1.0}]
y[0.5] /. sol
{1.23391}

To apply Taylor's algorithm, we pick up a fixed step length h = 0.1, which means that we need to do 10 iterations to reach the end point x = 1.

Clear[y]
h:=0.1;
f1[x_, y_] := 1 / (2 x - 3 y + 5)
fx[x_, y_] := D[f1[a, b], a] /. {a -> x, b -> y}
fy[x_, y_] := D[f1[a, b], b] /. {a -> x, b -> y}
f2[x_, y_] := fx[x, y] + fy[x, y] * f1[x, y]
y[0] = 1;
Do[y[n + 1] = y[n] + h * f1[h n, y[n]] + h^2 * f2[h n, y[n]]/2, {n, 0, 9}]
y[10]
Out[29]= 1.43529

At each step in the sequence, the expansion requires 3 addition, 1 subtraction, 3 multiplication, and 1 division operations. Since there are 8 steps to this sequence, in order to obtain the desired approximation, the code performs 80 operations, which would be hectic work to do by hand. The approximation is y10 = 1.435290196. This algorithm requires at least 22 operations per step, which means that the entire sequence requires 220 steps. However, this is not significant since Mathematica does all the work.

To check the answer, we write a subroutine to evaluate numerically the initial value problem for arbitrary slope function, initial conditions, and step size:

taylor2[f_, {x_, x0_, xn_}, {y_, y0_}, stepsize_] :=
Block[{xold = x0, yold = y0, sollist = {{x0, y0}}, h},
h = N[stepsize];
Do[xnew = xold + h;
fold = f /. {x -> xold, y -> yold};
Phi]2 = fold + (h/2)*((D[f, x] /. {x -> xold, y -> yold}) + (D[f, y] /. {x -> xold, y -> yold})*(f /. {x -> xold, y -> yold})); ynew = yold + h*\[Phi]2; sollist = Append[sollist, {xnew, ynew}]; xold = xnew; yold = ynew, {(xn - x0)/stepsize}]; Return[sollist]] Now we apply this subroutine to our slope function: taylor2[1/(2*x - 3*y + 5), {x, 0, 1}, {y, 1}, .1] // Flatten Out[8]= {0, 1, 0.1, 1.04938, 0.2, 1.09747, 0.3, 1.14427, 0.4, 1.18976, 0.5, 1.23393, 0.6, 1.27678, 0.7, 1.31833, 0.8, 1.35858, 0.9, 1.39756, 1., 1.43529} Using Mathematica command NDSolve, we check the obtained value: sol = NDSolve[{y'[x] == 1/(2*x - 3*y[x] + 5), y[0] == 1}, y[x], {x, 0, 2}] Out[5]= {{y[x] -> InterpolatingFunction[{{0., 2.}}<>][x]}} z = sol /. x -> sol /. x -> 1.0 Out[6]= {{y[{{y[1.] -> 1.43529}}] -> {{InterpolatingFunction[0.,2.}},<>][y[1.] -> 1.43529]}}}} The results of applying Taylor's approximation of order 2 with h = 0.1 n xn yn Exact solution Absolute error 1 0.1 1.04938 1.04937 -4.904×10-6 2 0.2 1.09747 1.09746 -8.9437×10-6 3 0.3 1.14427 1.14426 -0.000012028 4 0.4 1.18976 1.18974 -0.000014039 5 0.5 1.23393 1.23391 -0.0000149395 6 0.6 1.27678 1.27677 -0.0000147313 7 0.7 1.31833 1.31832 -0.0000134095 8 0.8 1.35858 1.35857 -0.0000110424 9 0.9 1.39756 1.39755 -7.70324×10-6 10 1.0 1.43529 1.43529 -3.49865×10-6 ■ Example: Let us take a look at the initial value problem $$y’=(x^2 -y^2) \sin (y), \quad y(0)=1 ,$$ for which it is impossible to find the actual solution. First, we apply the numerical solver to idenify the numerical values at grid points, sol = NDSolve[{Derivative[1][u][x] == (x^2 - u[x]^2)*Sin[u[x]], u[0] == 1}, u, {x, 0, 1.0}] u[0.9] /. sol Upon taking the mesh of points with uniform step size h = 0.1, $$x_n = n\,h , \quad n=0,1,2,\ldots ,$$ we denote by yn approximate values of unknown solution at these grid points. The initial point (0, 1) is obtained from the given initial condition. For second order Taylor's approximation, we need only to determine the first derivative of the slope function: \begin{align*} y' &= f(x,y) = \left( x^2 - y^2 \right) \sin (y) , \\ y'' &= 2\,x\,\sin y - 2\,y \left( x^2- y^2 \right) \sin^2 (y) + \left( x^2 - y^2 \right) \sin (y) \, \cos (y) . \end{align*} We check the formulas with Mathematica: f[x_, y_] = (x^2 - y^2)*Sin[y] f2[x_,y_] = D[f[x,y],x] + D[f[x,y],y] f[x,y] Then we ask Mathematica to do all calculations y[0] = 1; x[0] = 0; h = 0.1; Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2*f2[h*n, y[n]]/2, {n, 0, 10}] y[10] Here is the direct application of polynomial approximation with many loops. f[x_,y_]= (x^2 - y^2)*Sin[y] f1[x_,y_]=D[f[x,y],x]+f[x,y]*D[f[x,y],y] y[1] = 1; h =0.1; a = 0; b = 1; n = Floor[(b-a)/h]; Do[x[i] = 0.0 + (i-1)*h, {i,1,n+1}] Do[{k1=f1[x[i],y[i]], y[i+1] = y[i] + h*f[x[i],y[i]] +h*h/2*k1}, {i,1,n}] Do[Print[x[i], " ", y[i]],{i,1,n+1}] The results of applying Taylor's approximation of order 2 with h = 0.1 n xn yn Exact solution Absolute error 1 0.1 0.925207 0.924517 -0.000689672 2 0.2 0.865145 0.864125 -0.00101984 3 0.3 0.817529 0.816377 -0.0011522 4 0.4 0.779706 -0.0011734 5 0.5 0.754281 0.753152 -0.0011288 6 0.6 0.737247 0.736205 -0.00104159 7 0.7 0.729638 0.728715 -0.000922592 8 0.8 0.731621 0.730845 -0.000776388 9 0.9 0.743648 0.743043 -0.000605431 10 1.0 0.766427 0.76601293 -0.000414068 ■ Third Order Polynomial Approximations The third order Taylor approximation is obtained by adding a third order differential deviation to the equation for the second order expansion. \[ y_{n+1} = y_n + h\,f (x_n , y_n ) + \frac{h^2}{2} \ y'' (_n ) + \frac{h^3}{3!}\, y''' (x_n ) = y_n + h\, \Phi_3 (h) ,
where the increment function Φ3 adds just one term to Φ2:
$\Phi_3 (h) = f (x_n , y_n ) + \frac{h}{2} \left. \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right) f(x,y) \right\vert_{x=x_n , y= y_n} + \frac{h^2}{6} \left. \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right)^2 f(x,y) \right\vert_{x=x_n , y= y_n}$
Naturally, we expect the values defined by the above Taylor's algorithm to approximate the corresponding values y(xn) of the exact solution, particularly when step size h is small.

We will show how it works in two examples that were considered previously. As usual, we use a uniform step size h.

Example: We reconsider the initial value problem for the Riccati equation
$y' = y^2 - x^2 , \qquad y(0) = 1/2.$
Upon taking the mesh of points with step size h = 0.1, $$x_n = n\,h , \quad n=0,1,2,\ldots ,$$ we denote by yn approximate values of unknown solution at these grid points. The initial point (0, 0.5) is obtained from the given initial condition. In order to apply the Taylor approximation, we need to determine the derivatives of the slope function:
\begin{align*} f' (x,y) &= \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right) f(x,y) = 2y \left( y^2 - x^2 \right) - 2x = 2\left( y^3 - y\,x^2 - x \right) , \\ f'' (x,y) &= 6\,y^2 \left( y^2 - x^2 \right) -2\,x^2 \left( y^2 - x^2 \right) -4\,y\,x -2 = 6\,y^4 -8\,x^2 y^2 -2\,x^4 - 4\,y\,x -2 . \end{align*}
This tedious job is delegate to Mathematica:
f[x_,y_] = y^2 - x^2
op[x_,y_] = Expand[D[#,x] + f[x,y]*D[#,y]]&
Then we evaluate two derivatives:
phi1[x_,y_] = op[x,y][f[x,y]]
phi2[x_,y_] = op[x,y][op[x,y][f[x,y]]]
Now we are ready to perform calculations. For simplicity, we choose a constant step size
y[0] = 0.5; h = 0.1; x[0] = 0;
Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2/2*phi1[n*h, y[n]] + 1/6*h^3*phi2[h*n, y[n]] , {n, 0, 5}] y[5]
This gives the Taylor approximation of order three to be 0.617834, although the correct value is 0.6177691815444183. However, if you try step size h = 0.5 and make one step, the corresponding approximation becomes 0.622396. With smaller step size h = 0.01, we get 0.617769, but it requires 50 steps.

Now we calculate how many arithmetic operations are need to perform the third order Taylor approximation in this particular example. To evaluate f(x,y) requirs 3 arithmetic operations (2 multiplications an d1 subruction). To evaluate phi1, 7 arithmetic operations are needed, and for phi2 we need 13 arithmetic operations. Adding all these numbers and using 11 arithmetic operations for the algorithm, we come with 34 arithmetic operations at every step.

The results of applying Taylor's approximation of order 3 with h = 0.1
n xn yn Exact solution Absolute error
1 0.1 0.525979 0.525974 -5.57311×10-6
2 0.2 0.552752 0.552738 -0.0000139251
3 0.3 0.578443 0.578417 -0.0000259783
4 0.4 0.600943 0.600901 -0.0000426144
5 0.5 0.617834 0.617769 -0.0000644693
6 0.6 0.626319 0.626228 -0.0000915256
7 0.7 0.623174 0.623052 -0.000122428
8 0.8 0.604728 0.604575 -0.000153514
9 0.9 0.566941 0.566763 -0.000177769
10 1.0 0.505615 0.505431 -0.000184346
■
Example: We reconsider the IVP: $$y’=1/(2x-3y+5), \quad y(0)=1 .$$ To use the third order polynomial approximation, we choose a uniform grid with constant step length h = 0.1. In order to apply the Taylor approximation, we need to determine the derivatives of the slope function:
\begin{align*} f' (x,y) &= \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right) f(x,y) = \frac{3 - 2 \left( 2x-3y+5 \right)}{\left( 2x-3y+5 \right)^3} \\ f'' (x,y) &= \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right)^2 f(x,y) = \frac{27 -30\left( 2x-3y+5 \right)}{\left( 2x-3y+5 \right)^5} + \frac{8}{\left( 2x-3y+5 \right)^3} . \end{align*}

Here is the syntax for the third order Taylor approximation:

f[x_, y_] = 1/(2*x - 3*y + 5)
op[x_, y_] = Expand[D[#, x] + f[x, y]*D[#, y]] &
phi1[x_, y_] = op[x, y][f[x, y]]
phi2[x_, y_] = op[x, y][op[x, y][f[x, y]]]

y[0] = 1 ;h = 0.1; x[0] = 0
Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2*phi1[h*n, y[n]]/2 + h^3*phi2[h*n, y[n]]/6, {n, 0, 10}]
y[10]
Out[41]= 1.43529
The approximation for this value is y10=1.435283295 while the true value is 1.4352866048707.

The total number of numerical operations here is 420. It is obvious at this point why using mathematical programs such as Mathematica is a desired approach for such problems.

We can also apply the same approach from the previous example.
The results of applying Taylor's approximation of order 3 with h = 0.1
n xn yn Exact solution Absolute error
1 0.1 1.04938 1.04937 -4.904×10-6
2 0.2 1.09747 1.09746 -8.9437×10-6
3 0.3 1.14427 1.14426 -0.000012028
4 0.4 1.18976 1.18974 -0.0000140392
5 0.5 1.23393 1.23391 -0.0000149395
6 0.6 1.27678 1.27677 -0.0000147313
7 0.7 1.31833 1.31832 -0.0000134095
8 0.8 1.35858 1.35857 -0.0000110424
9 0.9 1.39756 1.39755 -7.70324×10-6
10 1.0 1.43529 1.43529 -3.49865×10-6
■

Example: We reconsider the IVP: $$y’=(x^2 -y^2) \sin (y), \quad y(0)=1 .$$ To use the third order polynomial approximation, we choose a uniform grid with constant step length h = 0.1.

taylor3[f_, {x_, x0_, xn_}, {y_, y0_}, stepsize_] :=
Block[{xold = x0, yold = y0, sollist = {{x0, y0}}, h},
h = N[stepsize];
Do[xnew = xold + h;
fold = f /. {x -> xold, y -> yold};
$Phi]2 = fold + (h/2)*((D[f, x] /. {x -> xold, y -> yold}) + (D[f, y] /. {x -> xold, y -> yold})*(f /. {x -> xold, y -> yold})); \[Alpha] = (D[f, {x, 2}] /. {x -> xold, y -> yold}) + 2*(D[f, x, y] /. {x -> xold, y -> yold})*(f /. {x -> xold, y -> yold}) + (D[f, {y, 2}] /. {x -> xold, y -> yold})*(f /. {x -> xold, y -> yold})^2 + (D[f, x] /. {x -> xold, y -> yold})*(D[f, y] /. {x -> xold, y -> yold}) + (D[f, y] /. {x -> xold, y -> yold})^2*(f /. {x -> xold, y -> yold}); ynew = yold + h*\[Phi]2 + h^3/3!*\[Alpha]; sollist = Append[sollist, {xnew, ynew}]; xold = xnew; yold = ynew, {(xn - x0)/stepsize}]; Return[sollist]] Solution = taylor3[(x^2 - y^2)*Sin[y], {x, 0, 1}, {y, 1}, .1] Out[8]= {{0, 1}, {0.1, 1.09483}, {0.2, 1.21862}, {0.3, 1.38695}, {0.4, 1.62303}, {0.5, 1.95295}, {0.6, 2.36595}, {0.7, 2.73957}, {0.8, 2.9658}, {0.9, 3.07578}, {1., 3.12003}} Another example of code (with the same output): f[x_,y_]=(x^2 + y^2)*Sin[y] f1[x_,y_]=D[f[x,y],x]+f[x,y]*D[f[x,y],y] f2[x_,y_]=D[f[x,y],x,x]+2*f[x,y]*D[f[x,y],x,y] + D[f[x,y],y,y]*f[x,y]*f[x,y] + D[f[x,y],x]*D[f[x,y],y] +f[x,y]*(D[f[x,y],y])^2 y[1] = 1; h =0.1; a = 0; b = 1; n = Floor[(b-a)/h]; Do[x[i] = 0.0 + (i-1)*h, {i,1,n+1}] Do[{k1=f1[x[i],y[i]], k2=f2[x[i],y[i]], y[i+1] = y[i] + h*f[x[i],y[i]] +h*h/2*k1 + k2/6*h^3}, {i,1,n}] Do[Print[x[i], " ", y[i]],{i,1,n+1}] The results of applying Taylor's approximation of order 3 with h = 0.1 n xn yn Exact solution Absolute error 1 0.1 0.92444 0.924517 0.0000773961 2 0.2 0.864011 0.864125 0.00011362 3 0.3 0.816248 0.816377 0.000128852 4 0.4 0.77957 0.779706 0.000133904 5 0.5 0.753018 0.753152 0.000134541 6 0.6 0.736072 0.736205 0.000133674 7 0.7 0.728582 0.728715 0.000132868 8 0.8 0.730712 0.730845 0.000132538 9 0.9 0.74291 0.743043 0.000132268 10 1.0 0.765882 0.76601293 0.000130676 ■ Fourth Order Polynomial Approximations To solve numerically the initial value problem \[ y' = f(x,y) , \qquad y(x_0 ) = y_0 ,$
we introduce the grid points starting with the initial point: x0, x1, … , xN, where, for simplicity, the step size is assumed to be uniform. Then we express these points according to the formula xn = x0 + n h, n = 0,1, … , N, with fixed step size h. We shall always denote by yn a number intended to approximate y(xn), the value of the true solution of the given initial value problem at x = xn. If the method of Taylor expansion of order 4 is used, these values are calculated according to the following scheme:
$y_{n+1} = y_n + h\,\Phi_4 (h) ,$
where
$\Phi_4 (h) = f(x,y) + \left. \frac{h}{2} \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right) f(x,y) \right\vert_{x=x_n, y=y_n} + \left. \frac{h^2}{3!} \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right)^2 f(x,y) \right\vert_{x=x_n, y=y_n} + \left. \frac{h^3}{4!} \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right)^3 f(x,y) \right\vert_{x=x_n, y=y_n} .$
This scheme is known as the Taylor algorithm of order 4. It approximates the solution of the initial value problem $$y' = f(x,y), \ y(a ) = y_0$$ over the interval [𝑎,b] by evaluating derivatives y'', y''', and y'''' and using the Taylor polynomial at each step of size h.
Example: We reconsider the initial value problem for the Riccati equation
$y' = y^2 - x^2 , \qquad y(0) = 1/2.$
Upon taking the mesh of points with step size h = 0.1, $$x_n = n\,h , \quad n=0,1,2,\ldots ,$$ we denote by yn approximate values of unknown solution at these grid points. The initial point is (0, 0.5) is obtained from the given initial condition. In order to apply the Taylor approximation, we need to determine the derivatives of the slope function f(x,y) = y² - x² up to the order 3:
\begin{align*} f' (x,y) &= \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right) f(x,y) = 2y \left( y^2 - x^2 \right) - 2x = 2\,y^3 - 2\,y\,x^2 - 2\,x , \\ f'' (x,y) &= 6\,y^2 \left( y^2 - x^2 \right) -2\,x^2 \left( y^2 - x^2 \right) -4\,y\,x -2 = 6\,y^4 -8\,x^2 y^2 -2\,x^4 - 4\,y\,x -2 , \\ f''' (x,y) &= \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right)^3 f(x,y) = 4 \left( 3\,x^3 -y + 4\,x^4 y -5\, x\,y^2 -10\,x^2 y^3 + 6\,y^5 \right) . \end{align*}
This tedious job we delegate to Mathematica:
f[x_,y_] = y^2 - x^2
op[x_,y_] = Expand[D[#,x] + f[x,y]*D[#,y]]&
Then we evaluate two derivatives:
phi1[x_,y_] = op[x,y][f[x,y]]
phi2[x_,y_] = op[x,y][op[x,y][f[x,y]]]]
phi3[x_, y_] = op[x, y][op[x, y][op[x, y][f[x, y]]]]
-2 x - 2 x^2 y + 2 y^3 -2 + 2 x^4 - 4 x y - 8 x^2 y^2 + 6 y^4 12 x^3 - 4 y + 16 x^4 y - 20 x y^2 - 40 x^2 y^3 + 24 y^5
Now we are ready to perform calculations. For simplicity, we choose a constant step size h.
y[0] = 0.5; h = 0.1; x[0] = 0;
Do[y[n + 1] = y[n] + h f[n*h, y[n]] + h^2/2*phi1[n*h, y[n]] + 1/6*h^3*phi2[h*n, y[n]]+ 1/24*h^4*phi3[h*n, y[n]], {n, 0, 5}]
Entering y[5], we get the fourth order Taylor approximation to be 0.617772; compare with the true value: 0.6177691815444183.     ■
Example: Consider the initial value problem $$y’=1/(2x-3y+5), \quad y(0)=1$$ and let us try to find out y(1). In order to apply the Taylor approximation of order four, we need to determine the derivatives of the slope function:
\begin{align*} f' (x,y) &= \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right) f(x,y) = \frac{3 - 2 \left( 2x-3y+5 \right)}{\left( 2x-3y+5 \right)^3} \\ f'' (x,y) &= \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right)^2 f(x,y) = \frac{27 -30\left( 2x-3y+5 \right)}{\left( 2x-3y+5 \right)^5} + \frac{8}{\left( 2x-3y+5 \right)^3} , \\ f''' (x,y) &= \left( \frac{\partial}{\partial x} + f(x,y)\,\frac{\partial}{\partial y} \right)^3 f(x,y) = \frac{405}{\left( 2x-3y+5 \right)^7} - \frac{630}{\left( 2x-3y+5 \right)^6} + \frac{312}{\left( 2x-3y+5 \right)^5} - \frac{48}{\left( 2x-3y+5 \right)^4} . \end{align*}

We set for simplicity a fixed step length h = 0.1, which means that we need to do 10 iterations to reach the end point x = 1.

f[x_, y_] = 1/(2*x - 3*y + 5))
op[x_, y_] = Expand[D[#, x] + f[x, y]*D[#, y]] &
phi1[x_, y_] = op[x, y][f[x, y]]
phi2[x_, y_] = op[x, y][op[x, y][f[x, y]]]]
phi3[x_, y_] = op[x, y][op[x, y][op[x, y][f[x, y]]]]]
]
y[0] = 1.0; h = 0.1; x[0] = 0;]
Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2*phi1[h*n, y[n]]/2 + h^3*phi2[h*n, y[n]]/6 + h^4*phi3[h*n, y[n]]/24 , {n, 0, 10}]]
y[10]
Out[29]= 1.43529

At each step in the sequence, the expansion requires 3 addition, 1 subtraction, 3 multiplication, and 1 division operations. Since there are 8 steps to this sequence, in order to obtain the desired approximation, the code performs 80 operations, which would be hectic work to do by hand. The approximation is y10 = 1.435290196. This algorithm requires at least 22 operations per step, which means that the entire sequence requires 220 steps. However, this is not significant since Mathematica does all the work.

■
Example: Let us take a look at the initial value problem $$y’=(x^2 -y^2) \sin (y), \quad y(0)=1 ,$$ for which it is impossible to find the actual solution. Hence, we use a standard Mathematica numerical solver
sol = NDSolve[{Derivative[1][u][x] == (x^2 - u[x]^2)*Sin[u[x]], u[0] == 1}, u, {x, 0, 1.0}] u[1.0] /. sol
0.766012930841635
First we find three derivatives:
\begin{align*} y' &= f(x,y) = \left( x^2 - y^2 \right) \sin (y) , \\ y'' &= 2\left( x- y\, y' \right) \sin (y) + \left( x^2 - y^2 \right) y' \, \cos (y) , \\ y''' &= \cdots , \end{align*}
where we did not type all 38 terms. This tedious job we delegate to Mathematica:
f[x_, y_] = (x^2 - y^2)*Sin[y]
op[x_, y_] = Expand[D[#, x] + f[x, y]*D[#, y]] & phi1[x_, y_] = op[x, y][f[x, y]] phi2[x_, y_] = op[x, y][op[x, y][f[x, y]]] phi3[x_, y_] = op[x, y][op[x, y][op[x, y][f[x, y]]]]
Then we ask Mathematica to do all calculations
y[0] = 1; x[0] = 0; h = 0.1;
Do[y[n + 1] = y[n] + h*f[n*h, y[n]] + h^2/2*phi1[n*h, y[n]] + 1/6*h^3*phi2[h*n, y[n]] + 1/24*h^4*phi3[h*n, y[n]], {n, 0, 10}]
Do[Print[y[n]], {n, 10}]
We compare the fourth order Taylor series approximation y10 = 0.766021 with the true value that gives an error 8.08804×10-6. So we see that Taylor's approximation of order 4 gives 5 correct significant digits.     ■