Numerical Approximation of Differential Equations with Midpoint Method I want to proof that the local truncation error of the Midpoint Method is  
$d_{k+1}=O\left(h^{3}\right)$

Approach
The local truncation error is defined as:
$d_{k+1}=u(t_{k+1})-\underset{u_{k+1}}{\underbrace{\left[u(t_{k})+\triangle t\cdot\phi\right]}}$


*

*$u(t_{k+1})\rightarrow$  exact solution of the differential equation $u'(t)=f(u,t)$, which is being approximated

*$\Phi\rightarrow f\left(u(t_{k+1})+\frac{1}{2}\left(\triangle t\right)\cdot f\left(u(t_{k}),t_{k}\right),t_{k}+\frac{1}{2}\left(\triangle t\right)\right)$ incremental Function called the Midpoint Method
The first step is to write the term $u(t_{k+1})$ as its Taylor expansion:
$u(t_{k+1})=u(t_{k})+\frac{1}{1!}\cdot h^{1}\cdot u^{(1)}(t_{k})+\frac{1}{2!}\cdot h^{2}\cdot u^{(2)}(t_{k})+\frac{1}{3!}\cdot h^{3}\cdot u^{(3)}(t_{k})+O\left(h^{4}\right)$
Then the second step is to write the incremental function in a fashion that the elements can cancel out with the previously introduced Taylor expansion of $u(t_{k+1})$.
I found a solution to the second step, but i do not understand how it is formed. And i would be happy if someone could explain in detail and write out the rules which 
lead to it. The solution goes as follows:
$\Phi\rightarrow f\left(u(t_{k+1})+\frac{1}{2}\left(\triangle t\right)\cdot f\left(u(t_{k}),t_{k}\right),t_{k}+\frac{1}{2}\left(\triangle t\right)\right)$
What happened here?
$=f\left(u\left(t_{k}\right),t_{k}\right)+\frac{1}{2}\cdot h\cdot f_{t}\left(u\left(t_{k}\right),t_{k}\right)+\frac{1}{2}\cdot h\cdot f\left(u\left(t_{k}\right),t_{k}\right)\cdot f_{y}\left(u\left(t_{k}\right),t_{k}\right)$
$+\frac{1}{2}\cdot\left(\frac{1}{2}\cdot h\right)^{2}\cdot f_{tt}+\left(\frac{1}{2}\cdot h\right)^{2}f\cdot f_{tu}+\frac{1}{2}\left(\frac{1}{2}\cdot h\right)^{2}\cdot f^{2}f_{uu}+O\left(h^{3}\right)$
I understand that there is some kind of differentiating going on, but i cant figure it out. Can someone provide the steps involved?

Edits
I found out that total derivatives are used somehow, but i do not have written out the exact transformation steps.
 A: Too much going on. Keep it simple to not get lost. Just look at the Taylor expansion about the midpoint. Then since one side is larger than the midpoint and one side is smaller, each with a distance of $\frac{\Delta t}{2}$, you get
$$ u(t+\Delta t) = u(t+\frac{\Delta t}{2})  + (\frac{\Delta t}{2})u'(t+ \frac{\Delta t}{2}) + (\frac{\Delta t}{2})^2 u''(t + \frac{\Delta t}{2}) + \mathcal{O}(\Delta t)^3
$$
Here I use shorthand $f(t)=f(u(t),t)$ and plugging in $u' = f$ you get 
$$u(t+\Delta t) = u(t+\frac{\Delta t}{2})  + (\frac{\Delta t}{2})f(t+ \frac{\Delta t}{2}) + (\frac{\Delta t}{2})^2 f'(t + \frac{\Delta t}{2}) +  \mathcal{O}(\Delta t)^3 $$
Call that equation 1. The same steps for the other side gives
$$u(t) = u(t+\frac{\Delta t}{2})  - (\frac{\Delta t}{2})f(t + \frac{\Delta t}{2}) + (\frac{\Delta t}{2})^2 f'(t + \frac{\Delta t}{2}) +  \mathcal{O}(\Delta t)^3 $$
Call that equation 2. This was the only hard part. Now subtract equation 1 from equation 2 to make equation 3:
$$u(t+\Delta t) - u(t) =  (\Delta t) f(t + \frac{\Delta t}{2}) + \mathcal{O}(\Delta t)^3 $$
What this tells us is that if you know the solution at the midpoint and use it to calculate the next step then you have an error of $\mathcal{O}(\Delta t)^3$ which is what we want. But how do we find out $u$ at the midpoint, i.e. $u(t+\frac{\Delta t}{2})$? Add equation 1 to equation 2 and divide by 2:
$$\frac{u(t+\Delta t) +u(t)}{2} = u(t+\frac{\Delta t}{2}) + \mathcal{O}(\Delta t)^2$$
Notice the $\mathcal{O}(\Delta t)$ term cancelled. Plugging this into $f$ as our approximation for $u$ at the midpoint, we get
$$ f(t+\frac{\Delta t}{2},u(t+\frac{\Delta t}{2})) = f(t+\frac{\Delta t}{2},\frac{u(t+\Delta t) +u(t)}{2}) + \mathcal{O}(\Delta t)^2 $$
Notice that this is an order 2 approximation to $f$ at the midpoint, but since $f$ is multiplied by $\Delta t$, we get that the error is order 3, that is by substitution into equation 3 we get:
$$u(t+\Delta t) - u(t) =  (\Delta t) f(t+\frac{\Delta t}{2},\frac{u(t+\Delta t) +u(t)}{2}) + \mathcal{O}(\Delta t)^3 $$
which is the midpoint method with the error. I hope this illuminates "why" it works: the Taylor expansion at the center makes the order 1 term flip signs and cancel. This trick is used a lot in numerical methods.
A: It is easier to do this for an autonomous DE and then only for the scalar case, the error terms for the other cases have a similar structure. Using $u'(t)=f(u(t))$ and thus $u''(t)=f'(u(t))u'(t)=f'(u(t))f(u(t))$ etc. gives
\begin{align}
u(t+h)&=u(t)+u'(t)·h+\frac12·u''(t)·h^2+\frac16·u'''(t)·h^3+\frac1{24}·u^{(4)}(t)·h^4+…
\\
&=u(t)+f·h+\frac12·f'·f·h^2+\frac16·(f''·f^2+f'^2·f)·h^3
\\
&\qquad\qquad+\frac1{24}·(f'''·f^3+4·f''f'f^2+f'^3·f)·h^4+…
\end{align}
and for the improved Euler or midpoint method at $u=u(t)$
\begin{align}
u_+&=u+h·f\left(u+\frac h2·f\right)
\\
&=u+h·\left(f+f'·\frac h2·f+\frac12·f''·\left(\frac h2·f\right)^2+\frac16·f'''·\left(\frac h2·f\right)^3\right)+…
\\
&=u+f·h+\frac12·f'·f·h^2+\frac18·f''·f^2·h^3+\frac1{48}·f'''·f^3·h^4+…
\end{align}
where we see several differences in the coefficient of $h^3$. 

See also the (Heun-)Ralston method. https://en.wikipedia.org/wiki/Runge%E2%80%93Kutta_methods#Usage

For a comprehensive methodology to treat the vector case, see Butcher's rooted trees as presented in the introduction here

In the vector case, the terms in the expansion of the midpoint step all have the form $f^{(k)}(u)[f(u),…,f(u)]_{k\text{ repetitions}}$. The terms in the Taylor series get more structure as the derivatives now are tensors, $f'$ is a linear operator, $f''$ a vector valued symmetric 2-form, $f'''$ a vector-valued symmetric 3-form etc.
\begin{align}
u'&=f(u)\\
u''&=f'(u)[f(u)]\\
u'''&=f''(u)[f(u),f(u)]+f'(u)[f'(u)[f(u)]]\\
u^{(4)}&=f'''(u)[f(u),f(u),f(u)]+3·f''(u)[f'(u)[f(u)],f(u)]+f'(u)[f''(u)[f(u),f(u)]]+f'(u)[f'(u)[f'(u)[f(u)]]]
\end{align}
From this one can see why rooted trees are a good idea to describe the single terms.

As to the last point of the question, this is a Taylor expansion that IMO hides more of the structure of the problem than helps to expose it. The structure is
$$
u_+=u+h·f(u+h·a,t+h·b),\text{ where } a=\frac12·f(u,t)\text{ and }b=\frac12
$$
For the error analysis, the first $f$ gets replaced by its Taylor expansion in $(u,t)$. Note that here $a,b$ are the variables, $u,t,f(u,t)$ are constants. Thus
$$
u_+=u+h·\Biggl(f(u,t)+h·\Bigl(f_u(u,t)·a+f_t(u,t)·b\Bigr)\\
+\frac{h^2}{2}·\Bigl(f_{uu}(u,t)·a^2+2f_{ut}(u,t)·ab+f_{tt}(u,t)·b^2\Bigr)\\
+\frac{h^3}6·\Bigl( f_{uuu}(u,t)·a^3+3f_{uut}(u,t)·a^2b+3f_{utt}(u,t)·ab^2+f_{ttt}(u,t)·b^3\Bigr) +…\Biggr)
$$
Now insert the terms for $a$ and $b$ to get the full form of the expansion.
A: Another variant (sufficiently unrelated to my previous answer) to approach the error goes first over an integral and does not rely as dominantly on Taylor formulas.
For an autonomous ODE $y'=f(y)$ with some exact solution with  $y(0)=y_0$, $f_0=f(y(0))$ one can also decompose the local error  using the Picard integral formulation of the ODE as follows:
\begin{align}
y(h)-y_0-hf(y_0+\tfrac h2 f_0)
&=\int_0^hf(y(s))ds - hf(y_0+\tfrac h2 f_0)
\\
&=\int_0^{h/2}\Bigl[(f\circ y)(\tfrac h2+s)-2(f∘y)(\tfrac h2)+(f∘y)(\tfrac h2-s)\Bigr]\,ds
\\
&\quad+h\Bigl[f\bigl(y(\tfrac h2)\bigr)-f\bigl(y_0+\tfrac h2 f_0\bigr)\Bigr]
\end{align}
The integrand of the first term is a second order difference for $f∘y$ which can be bounded by $C·s^2$ where $C$ is a bound for the second derivative $(f∘y)''=f''f^2+f'^2f$ which gives an $O(h^3)$ contribution of the integral from $s=0$ to $\frac h2$. The second term is proportional to the size of $y(\frac h2)-y_0-\frac h2f_0$ which as the error of an Euler step is $O(h^2)$ which with the leading factor results also in a contribution of $O(h^3)$.
