How to derive the simplified Newton iteration in the TR-BDF2 ODE integration scheme

The Problem

The TR-BDF2 explained in this paper [1], is quite a popular numerical scheme used to integrate $dot{y} = f(t,y)$, consistent of the following two stages:

begin{align}
y_{n+gamma} & = y_n + gamma frac{h}{2}left( f_n + f_{n+gamma} right) tag{1}
%
y_{n+1} & = frac{1}{gamma(2-gamma)}y_{n+gamma} – frac{(1-gamma)^2}{gamma(2-gamma)}y_n +
frac{1-gamma}{2-gamma}hf_{n+1} tag{2}
end{align}

The above are two implicit algebraic equations that need to be solved in each step. This is accomplished by performing a simplified Newton (chord) iteration, where in accordance to Hosea & Shampine’s paper [1] is performed by firstly letting $z=hf(t,y)$. As a result by substituting this into equations $(1)$ and $(2)$ above for the $k^{text{th}}$ Newton iteration, we have:

begin{align}
y_{n+gamma}^k & = left(y_n + frac{gamma}{2}z_nright) + frac{ gamma}
{2}z_{n+gamma}^k tag{3}
%
y_{n+1}^k & = left(y_n +frac{sqrt{2}}{4}z_n + frac{sqrt{2}}{4}z_{n+gamma} right) + frac{gamma}{2}z_{n+1}^k tag{4}
end{align}

From the above, the paper goes straight away to say that the Newton step and Newton iteration are as follows:

begin{align}
left(I-hfrac{ gamma}{2}frac{partial f}{partial y}right)Delta^k = hf(t_{n+gamma},y^k_{n+gamma}) – z_{n+gamma}^k,
text{and} z_{n+gamma}^{k+1} = z_{n+gamma}^k + Delta^k text{for equation (3)} tag{5}
%
left(I-hfrac{ gamma}{2}frac{partial f}{partial y}right)Delta^k = hf(t_{n+1},y^k_{n+1}) – z_{n+1}^k,
text{and} z_{n+1}^{k+1} = z_{n+1}^k + Delta^k text{for equation (4)} tag{6}
end{align}

The Question

Therefore I do not understand how the last two equations above were derived, and what the quantity $z$ represents.

To elaborate, the paper says that $z=hf(t,y)$, so $z_{n+gamma}^k=hf(t_{n+gamma},y_{n+gamma}^k)$, but doesn’t that make the RHS of equations $(5)$ and $(6)$ zero? Also, the typical Newton iteration is of the form $u^{k+1} = u^k – (g’)^{-1}g(u^k)$. In this case then what is the function $g$ corresponding to for equations $(5)$ and $(6)$? It surely can’t be $g = hf(t_{n+gamma},y^k_{n+gamma}) – z_{n+gamma}^k$, since the Jacobian of this function is not $I-hfrac{ gamma}{2}frac{partial f}{partial y}$ as shown in equation $(5)$ for example.

Thanks for your help in advance.

Bibliography

1. Hosea, M. E.; Shampine, L. F., Analysis and implementation of TR-BDF2, Appl. Numer. Math. 20, No. 1-2, 21-37 (1996). ZBL0859.65076.