El método de Euler se basa en la proximidad de la recta tangente a los valores de la función en un entorno del punto de tangencia, es por eso una aproximación lineal. Para aproximaciones de órdenes superiores, no lineales, puede recurrirse a los polinomios de Taylor.

El polinomios de Taylor de grado $p$ centrado en $x=x_0$ tiene la forma $$P_p(x)=y(x_0)+y'(x_0)(x-x_0)+\frac{y''(x_0)}{2!}(x-x_0)^2+\ldots+ \frac{y^{(p}(x_0)}{p!}(x-x_0)^p$$ Para escribir el polinomio de Taylor de la solución $y=y(x)$ del problema de valor inicial $$\left\{\begin{array}{ll} y'=f(x,y)\\ y(x_0)=y_0\end{array}\right.$$ se necesitarán las derivadas de $y(x)$ en $x_0$ desde orden 0 a orden $p$ :

• $y(x_0)$ viene dado por el problema

• $y'(x_0)$ se obtiene sustituyendo $x$ por $x_0$ en $y'(x)=f(x,y(x))$

• para hallar $y''(x_0)$ se deriva implícitamente $y'(x)=f(x,y(x))$ y se sustituye después $x$ por $x_0$; habrá que utilizar el valor hallado antes para $y'(x_0)$

• de la misma manera que esta última se irán calculando las siguientes, utilizando siempre los valores calculados en los pasos anteriores.

Si llamamos $h=x-x_0$, el polinomio anterior proporciona un método recursivo de aproximación de la solución $y(x)$ en un conjunto de puntos cercanos a $x_0$, que viene dada por las fórmulas

\begin{eqnarray*} x_{n+1}&=&x_n+h \\ y_{n+1}&=&y_n+hf(x_n,y_n)^+\frac{h^2}{2!}f_2(x_n,y_n)+\ldots+\frac{h^p}{p!}f_p(x_n,y_n)\end{eqnarray*} donde $f_k$ denota la expresión obtenida para la derivada $y^{(k}$ a partir de $y'=f(x,y)$.

Éste es el método de Taylor de orden $p$, del que se demuestra que tiene una velocidad de convergencia igual a $p$.