The equivalent integral equation for the time step is
$$
y(t+h)=y(t)+\int_0^hf(t+s,y(t+s))ds\tag{I}
$$
The overall method uses the trapezoidal formula
$$
\int_a^b u(s)ds=\frac{u(a)+u(b)}2(b-a)+O((b-a)^3)
$$
to find
$$
y_{n+1}=y_n+\frac h2(f(t_n,y_n)+f(t_{n+1},\hat y_{n+1})).
$$
Now the intermediate value $\hat y_{n+1}$ can be obtained in any way that is consistent, that is first order accurate, to get a second order method. Choosing the explicit Euler step gives the Heun method, setting $\hat y_{n+1}= y_{n+1}$ gives the implicit trapezoidal method, or to make it extra strange one can also select the value of the implicit Euler step as done here.
The implicit Euler step can as well be motivated by the integral formulation of the ODE, approximate (I) by the right-sided Riemann sum
$$
y(t+h)\approx y(t)+f(t+h,y(t+h))h
$$
to get the formula for the implicit Euler method.
What you have now in the composite method is a second order method similar to the trapezoidal method, the third order conditions fail already in the quadrature condition $b_1c_1^2+b_2c_2^2=\frac13$, as the left side is $\frac12$.
For the linear system $z'(t)=A(t)z(t)$ you can solve the stages directly, $k_1=A(t)z(t)$, $k_2=(I-A(t+h)h)^{-1}A(t+h)z(t)$. It is usually better to avoid the computation of the inverse matrix and just use a solver for linear systems.
