Your questions contains multiple subquestions.
- The advantage of the explicit method is that it requires only two function evaluations pr. time step. The advantage of the implicit method is that it is A-stable.
- The disadvantage of the explicit method is that the stability region is small. The disadvantage of the implicit method is that typically requires more than two function evaluations per time step.
- The stability of these method are studied by examining their behavior when applied to the simple test equation $y' = \lambda y$. When the real part of $\lambda$ is strictly less than zero, then the exact solution decays to zero as $t$ tends to infinity. The implicit method will reproduce this behavior regardless of the time step $h$. The explicit method will only produce this behavior for sufficient small values of $h$.
Using the implicit method requires less knowledge about the problem, but it will frequently take more time to obtain a solution. You can often reduce the time to solve the nonlinear equation by using the explicit method to generate a good initial guess. If you are dealing with one-dimensional problem, then the secant method can further reduce the solve time. If you have a multi-dimensional problem, then you should try to apply Newton's method.
When judging the quality of a scheme we should be concerned with the extent to which it reproduces the key features of the physical reality. There are several kinds of stability and which one is relevant to you depends on what kind of problem you are trying to solve.
Most textbooks which discuss stiff equations will have much more information about the terms "stability", "stability region", "A-stable methods". The simple test equation may appear trivial, but appearances are deceptive and it actually covers an amazing large class of physical problems.
Both Explicit Trapezoid method and Heun's method are 2nd order Runge Kutta methods in this form:
$$y_{n+1}=y_n+b_1hf(t_n,y_n)+b_2hf(t_i+c_2h,y_n+a_{21}hf(t_n,y_n)).$$
For them to be second order, they need to satisfy:
$$b_1+b_2=1, b_2c_1=1/2, a_{21}b_2=1/2.$$
For Explicit Trapezoid, we have
$$b_1=b_2=1/2, c_2=1, a_{21}=1.$$
For Heun's, we have
$$b_1=1/4, b_2=3/4, c_2=2/3, a_{21}=1.$$
For your example, Heun's method gives:
$$y_n+\frac{h}{4}\cdot\frac{8\cos(4t_n)+y_n}{4}+\frac{3h}{4}\frac{8\cos(4(t_n+\frac{2}{3}h))+y_n+h\cdot\frac{8\cos(4t_n)+y_n}{4}}{4}.$$
I think you can find this in any Numerical Analysis book.
Your Trapezoid method formula is correct, but you didn't apply it correctly in your example. It should be:
$$y_n+\frac{h}{2}\frac{8\cos(4t_n)+y_n}{4}+\frac{h}{2}\frac{8\cos(4(t_n+h))+h\frac{8\cos(4t_n)+y_n}{4}}{4}.$$
The 2nd order Taylor method is a totally different method. It uses Taylor expansion, so requires more partial derivatives of $f$:
$$y_{n+1}=y_n+hf(t_n,y_n)+\frac{h^2}{2}\left(\frac{\partial f}{\partial t}(t_n,y_n)+\frac{\partial f}{\partial y}(t_n,y_n)f(t_n,y_n)\right).$$
Applying to your example, we get
$$y_{n+1}=y_n+h\frac{8\cos(4t_n)+y_n}{4}+\frac{h^2}{2}\left(-8\sin(4t_n)+\frac{1}{4}\frac{8\cos(4t_n)+y_n}{4}\right).$$
Best Answer
The resulting problem would be: $$\left\{\begin{align*} \dot{x}&=u\\ \dot{y}&=v\\ \dot{u}&=-x/r^3\\ \dot{v}&=-y/r^3 \end{align*}\right.$$ with $u(0) = v(0) =y(0)= 0$ and $x(0) = 0.5$
Now you have with $\vec{w} = (x, y, u, v)^{T}$ and $\vec{f}=(u,v,-x/r^3,-y/r^3)^{T}$the following system: $$\frac{d\vec{w}}{dt} = \vec{f}(\vec{w})\tag{*}$$ Simply apply your integration rule to $(*)$. Since your numerical method is implicit, and $\vec{f}$ is non linear in $\vec{w}$, you will need some iterative solver to obtain the $\vec{w}^{n+1}$ vector solution.
Going further, permit me to apply the method to $(*)$. The resulting numerical scheme is: $$\vec{w}^{n+1}=\vec{w}^n+\frac{h}{2}\left(\vec{f}(\vec{w}^{n})+\vec{f}(\vec{w}^{n+1})\right)$$ The term $\vec{f}$ this time does not depend explicitly on time $t$.