I'm reading J. Betts' book "Practical Methods for Optimal Control and Estimation Using Nonlinear Programming" and in the chapter on IVP there is Butcher tableau for Hermite-Simpson method: $$ \begin{array}{c|cccc} 0 & 0 & 0 & 0 \\ \frac{1}{2} & \frac{5}{24} & \frac{1}{3} & -\frac{1}{24} \\ 1 & \frac{1}{6} & \frac{2}{3} & \frac{1}{6} \\ \hline & \frac{1}{6} & \frac{2}{3} & \frac{1}{6}\end{array} $$
General notation used is the following: $$ \begin{array}{c|cccc} \rho_1 & \alpha_{11} & \ldots & \alpha_{1K} \\ \vdots & \vdots & & \vdots \\ \rho_K & \alpha_{K1} & \ldots & \alpha_{KK} \\ \hline & \beta_1 & \ldots & \beta_k \end{array} $$
Also there are formulas: $$ \begin{align*} &\mathbf{y}_{i+1} = \mathbf{y}_i + h_i\sum\limits_{j=1}^K \beta_j \mathbf{f}_{ij} \\ & \text{where } \mathbf{f}_{ij} = \mathbf{f}\left[\left(\mathbf{y}_i + h_i\sum\limits_{l=1}^K \alpha_{jl}\mathbf{f}_{il} \right), (t_i + h_i\rho_j)\right] \end{align*} $$
I don't understand how to get formula for intermediate point $\mathbf{y}_{i+1/2} = \frac{1}{2}(\mathbf{y}_i + \mathbf{y}_{i+1}) + \frac{h_i}{8}(\mathbf{f}_i - \mathbf{f}_{i+1})$ using these formulas above.
So I use the coefficients from the second row and get some ugly formula: $$ \mathbf{y}_{i+1/2} = \mathbf{y}_{i} + h_i \left(\frac{5}{24}\mathbf{f}_{i1} + \frac{1}{3}\mathbf{f}_{i2} - \frac{1}{24}\mathbf{f}_{i3}\right) $$ and I completely stuck, the only thing I managed to see is that $\mathbf{f}_{i1} = \mathbf{f}(\mathbf{y}_i, t_i) = \mathbf{f}_i$. But it doesn't help much. I need help because I'm confused
Using a little different notation, let's rewrite it in the following form: \begin{align} \mathbf{y}_{i+1/2} &= \mathbf{y}_{i} + h_i\left(\frac{5}{24}\mathbf{f}_{i} + \frac{1}{3}\mathbf{f}_{i+1/2} - \frac{1}{24}\mathbf{f}_{i+1}\right) \\ &= \mathbf{y}_{i} + \frac{h_i}{2}\left(\frac{5}{12}\mathbf{f}_{i} + \frac{2}{3}\mathbf{f}_{i+1/2} - \frac{1}{12}\mathbf{f}_{i+1}\right) \end{align}
Let's now notice that $\frac{5}{12}\mathbf{f}_{i} = \frac{1}{6}\mathbf{f}_{i} + \frac{1}{4}\mathbf{f}_{i}$ and $-\frac{1}{12}\mathbf{f}_{i+1}=\frac{1}{6}\mathbf{f}_{i+1} - \frac{1}{4}\mathbf{f}_{i+1}$. So we get: $$ \mathbf{y}_{i+1/2} = \mathbf{y}_{i} + \frac{h_i}{2}\left(\frac{1}{6}\mathbf{f}_{i} + \frac{2}{3}\mathbf{f}_{i+1/2} + \frac{1}{6}\mathbf{f}_{i+1}\right) + \frac{h_i}{8}\left(\mathbf{f}_i + \mathbf{f}_{i+1}\right) $$
Now let's notice that: $$ \mathbf{y}_{i+1/2} = \frac{\mathbf{y_i}}{2} + \frac{1}{2}\left(\underbrace{\mathbf{y}_i + h_i\left(\frac{1}{6}\mathbf{f}_{i} + \frac{2}{3}\mathbf{f}_{i+1/2} + \frac{1}{6}\mathbf{f}_{i+1}\right)}_{\mathbf{y}_{i+1}}\right) + \frac{h_i}{8}\left(\mathbf{f}_i + \mathbf{f}_{i+1}\right) $$
And finaly $$ \mathbf{y}_{i+1/2} = \frac{1}{2}\left(\mathbf{y}_{i} + \mathbf{y}_{i+1}\right) + \frac{h_i}{8}\left(\mathbf{f}_i + \mathbf{f}_{i+1}\right) $$