In my derivation for $y^\prime = f(t,y)$, I begin by writing the method as an expression which should simplify to the error, by substitution of the exact solution \begin{equation} y(t_{n+1}) - y(t_n) - hf\left(t_n+\frac{1}{2}h,\,\frac{1}{2}(y(t_{n+1}) + y(t_n))\right)\end{equation} From here, it may be simplified in a straightforward way using Taylor expansions to the following \begin{equation}hy^\prime(t_n) + \frac{1}{2}h^2y^{\prime\prime}(t_n) + \mathcal{O}(h^3) - hf\left(t_n+\frac{1}{2}h,\,y\left(t_{n}+\frac{1}{2}h\right)+ \mathcal{O}(h^2)\right)\end{equation} This is where I am unsure how to proceed. From this question, the top answer seems to suggest that it is a trivial matter to assume that we may assert $$f(t_n + 0.5h,y(t_n + 0.5h) + \mathcal{O}(h^2)) = f(t_n + 0.5h,y(t_n + 0.5h)) + \mathcal{O}(h^2) = y^\prime(t_n+0.5h) + \mathcal{O}(h^2)$$ but this does not seem obvious to me. While it is stipulated that $f$ is an analytic function, I don't see how a perturbation in the input would necessarily perturb the output on the order of $h^2$ in general.
If this is a good assumption to make, could someone please explain? If not, why? What is the correct way to proceed here?
In general you assume the Lipschitz condition, from there you get $|f(t,y+O(h^2))-f(t,y)|=O(h^2)$.
For an order $p$ method to get the order $p$ local error you also need in general that $f$ is $p$ times continuously differentiable in all variables. Or at least $p-1$ times differentiable, with the highest derivative Lipschitz continuous.
You should probably also use the central difference quotient for the first terms, $$ y(t+h)-y(t)=hy'(t+0.5h)+O(h^3), $$ then there are no second order terms in the formula that have to be compensated for.