 Solution

The nodes are $c=(0,1/2,1{)}^T$. Their Lagrange interpolation basis on the normalized step isIntegrating these polynomials from zero to each node gives $a_{ij}={\int }_0^{c_i}{\ell }_j(s)ds$. At $c_1=0$ the row vanishes; at $c_2=1/2$ the integrals are $(5/24,1/3,-1/24)$; at $c_3=1$ they are $(1/6,2/3,1/6)$. The endpoint integrals also give $b=(1/6,2/3,1/6{)}^T$. These are exactly the coefficients in the original PDF.

Indeed, the degree-three polynomial $q$ determined by $q(0)=y_n$ and $q^{\prime }(hs)={\sum }_j{\ell }_j(s)f(q(h{c}_j))$ satisfies $q^{\prime }(h{c}_i)=f(q(h{c}_i))$. Integrating it to the nodes gives precisely the stage equations, and integrating to the endpoint gives the update. This proves it is a collocation Runge-Kutta method, specifically the three-stage Lobatto IIIA method.

For an explicit order proof, the collocation integrals imply $A1=c$, $Ac=c^2/2$ and $A{c}^2=c^3/3$, with powers componentwise. The endpoint weights are Simpson weights, so $b^T{c}^j=1/(j+1)$ for $j=0,1,2,3$. Therefore the eight Butcher order conditions through order four are satisfied:For example the last expression is $b^T{c}^3/6$, because $A^{2}c=A{c}^2/2=c^3/6$. These conditions are obtained by matching all elementary differential terms through degree four in the exact and numerical Taylor expansions, not just by testing a scalar linear equation. Nonautonomous equations are included by adjoining $t^{\prime }=1$.

To exclude higher order, the scalar stability function computed below satisfies $R(z)-e^z=-z^5/720+O(z^6)$. The local defect is therefore nonzero at degree five even on the scalar test equation. Hence $\boxed{\text{the classical order is exactly }4}$.