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).$$
Basically the idea is that a general non-linear ODE can be linearized locally, given some reasonable constraints on the DE, and
Consider $y'(x) = f(y,x)$ with $y(x_0)=b$. Then consider:
$$ y'(x) \approx f(b,x)+\frac{\partial f}{\partial y}\bigg\rvert_{y=b} [y-b]\approx \lambda y(x) + f(b,x) - \lambda b = \lambda y + c $$
where we have linearized about $y=b$, fixed some $x=x_0$, and let $\lambda=\partial_y f(b,x_0)$.
Then of course we can solve the homogeneous ODE $y' =\lambda y$ to get a general solution locally.
The idea then is that if the algorithm has pleasant properties on this linearized version, then those properties carry over to the full ODE (as long as the ODE behaves nicely).
Here are some relevant links that go into more detail:
Best Answer
Stability of explicit Runge-Kutta methods around a fixed point or a slow moving equilibrium is obtained when $Lh\sim 1$ with $L$ the $y$-Lipschitz constant and $h$ the step size. The exact constant on the right side depends on the order of the method, usually it is between 1 and 5. For first and second order methods that is only a guideline as the stability region does not contain a semi-circle around the origin in the negative complex half-plane.
To get a practical feeling for this, just solve your problem for various step sizes over the integration interval from $0$ to $1/3$, as $y(1/3)=0.1e^{-50}$ is a really small number, relative to the initial value and floating point precision.
with the results
As the stability region suggests, $0<h<\frac2{150}$ are valid choices in the table above, giving results close to zero. Sensible results that are practially zero as floating point numbers (and relative to the initial value) are only obtained for about $h<\frac{1.5}{150}$.