Estimating the truncation error Consider the following method 
$y_{n+1}=y_n+hf(t_{n+1},y_n+hf(t_n,y_n))$, the truncation error is given by 
$T_n=\dfrac{y(t_{n+1})-y(t_n)}{h}-f(t_{n+1},y_n+hf(t_n,y_n)) $, Im having problem on how to estimate $T_n$ especially on how to apply Taylor series on $f(t_{n+1},y_n+hf(t_n,y_n))$. Any help I would appreciate.
suppose your $f$ to be an analytic function, then, given a numerical method $\mathcal{Y}(f,k,y(t_n))$, where $f$ is your function, $k$ is the step size (you call it $h$), and $y(t_n)$ is the analytical solution at the time $t_n$, then the local error is given by $y(t_{n+1})-\mathcal{Y}(f,k,y(t_n))$. In genereal therms, it's just the difference between the analytical solution at time $t_{n+1}$ and the numerical method applied to the analytical solution at time $t_n$.

Evaluating this difference, and remembering that $\vec f(t_n,y(t_n))=\vec y'(t_n)$, since we're considering first order systems of ODE, you have to center the Taylor series in the point $(t_n,y(t_n))$:
\begin{align}
&amp;y(t_{n+1})-\mathcal{Y}(f,k,y(t_n))
\\
&amp;= y(t_n) +ky'(t_n)+\frac{k^2}{2}y''(t_n) + \mathcal{O(k^3)}
\\&amp;\qquad-y(t_n)-k[y'(t_n)+k\frac{df}{dt}(t_{n+1},y(t_n)) + ky'(t_n)\frac{df}{dy}(t_{n},y(t_n))+\mathcal{O(k^2)}]
\\
&amp;=ky'(t_n)+\frac{k^2}{2}y''(t_n)+\mathcal{O(k^3)}-ky'(t_n)-k^2 \frac{df}{dt}(t_n,y(t_n))-k^2[y'(t_n) \frac{df}{dy}(t_n,y(t_n))]
\\
&amp;=ky'(t_n)+\frac{k^2}{2}y''(t_n)+\mathcal{O(k^3)}-ky'(t_n)-k^2 \frac{df}{dt}(t_n,y(t_n))-k^2[y''(t_n)-\frac{df}{dt}(t_n,y(t_n))] 
\\
&amp;\approx \mathcal{O(k^3)}
\end{align}
Since we're integrating in an interval of length $k$, then this method is of the second order. (And btw it's a particular Runge-Kutta method)

Hope it will be useful for you
suppose your $f$ to be an analytic function, then, given a numerical method $\mathcal{Y}(f,k,y(t_n))$, where $f$ is your function, $k$ is the step size (you call it $h$), and $y(t_n)$ is the analytical solution at the time $t_n$, then the local error is given by $y(t_{n+1})-\mathcal{Y}(f,k,y(t_n))$. In genereal therms, it's just the difference between the analytical solution at time $t_{n+1}$ and the numerical method applied to the analytical solution at time $t_n$.

Evaluating this difference, and remembering that $\vec f(t_n,y(t_n))=\vec y'(t_n)$, since we're considering first order systems of ODE, you have to center the Taylor series in the point $(t_n,y(t_n))$:
\begin{align}
&amp;y(t_{n+1})-\mathcal{Y}(f,k,y(t_n))
\\
&amp;= y(t_n) +ky'(t_n)+\frac{k^2}{2}y''(t_n) + \mathcal{O(k^3)}
\\&amp;\qquad-y(t_n)-k[y'(t_n)+k\frac{df}{dt}(t_{n+1},y(t_n)) + ky'(t_n)\frac{df}{dy}(t_{n},y(t_n))+\mathcal{O(k^2)}]
\\
&amp;=ky'(t_n)+\frac{k^2}{2}y''(t_n)+\mathcal{O(k^3)}-ky'(t_n)-k^2 \frac{df}{dt}(t_n,y(t_n))-k^2[y'(t_n) \frac{df}{dy}(t_n,y(t_n))]
\\
&amp;=ky'(t_n)+\frac{k^2}{2}y''(t_n)+\mathcal{O(k^3)}-ky'(t_n)-k^2 \frac{df}{dt}(t_n,y(t_n))-k^2[y''(t_n)-\frac{df}{dt}(t_n,y(t_n))] 
\\
&amp;\approx \mathcal{O(k^3)}
\end{align}
Since we're integrating in an interval of length $k$, then this method is of the second order. (And btw it's a particular Runge-Kutta method)

Hope it will be useful for you