There is a theorem which states the following.
Let $u:I⊂ℝ→ℝ$ be given by
\begin{align*}u^{(n)} &= F(t,u,u',…,u^{(n-1)})\\
u^{(k)}(t_0)&=u_k^0, \qquad k=0,…,n-1.
\end{align*}
This is called an explicit, scalar initial value problem (IVP) of order $n$. Then $u$ is also given as solution to the system of first order:
\begin{align*}
u'_{n-1}(t)&=F(t,u_0(t),u_1(t),…,u_{n-1}(t)), & u_{n-1}(t_0)&=u_{n-1}^0 \\
u'_{n-2}(t)&=u_{n-1}(t), & u_{n-2}(t_0)&=u_{n-2}^0 \\
\vdots& \\
u'_{0}(t)&=u_{1}(t),& u_{0}(t_0)&=u_{0}^0
\end{align*}
using the notation $u_0(t)=u(t)$ and $u_k=\frac{d^k}{dt^k}u$.
What this theorem states, is that you can rewrite a IVP of higher order to a system of ODEs of order one. Written as a matrix system, if $F$ is linear, you will get
$$u' = Au + b.$$
Note, that your ODE is in fact a linear ODE.
The reason why this theorem is important, is that all theories that are developed for order one IVPs, $u'=f(t,u)$, can be applied to such systems of order one. Especially you can use all numerical solvers like the Explicit Euler or RK-methods.
In your case (using your first result $(*)$) that equation is:
$$\begin{bmatrix} u_1' \\ u_2'\\u_3' \end{bmatrix} = \begin{bmatrix}0 & 1 & 0\\ 0& 0 & 1 \\ -1/x & -x &-1 \end{bmatrix} \begin{bmatrix}u_1 \\ u_2 \\ u_3 \end{bmatrix} +\begin{bmatrix}0\\0\\x\end{bmatrix}$$
Note, that you have to also include the equation like $u_1' = u_2$ in that system. Also note, that you can't have any $y$'s in there anymore.
That system could be solved by any ODE solver, e.g. Explicit Euler:
$$u^{(t+1)} = u^{(t)} + h(Au^{(t)}+b).$$
Is it clear how to set up your iteration from here?
Now for your other case $(**)$. You will get the following system, modulo my calculation errors:
$$\begin{bmatrix}1 & 0 & 0\\ 0& 1 & 0 \\ x^2 & x & x \end{bmatrix}\begin{bmatrix} u_1' \\ u_2'\\u_3' \end{bmatrix} = \begin{bmatrix}0 & 1 & 0\\ 0& 0 & 1 \\ -1 & 0 &0 \end{bmatrix} \begin{bmatrix}u_1 \\ u_2 \\ u_3 \end{bmatrix} +\begin{bmatrix}0\\0\\1\end{bmatrix}$$
So this is something like $$Mu'=\tilde{A}u+\tilde{b}.$$ If you would multiply by $M^{-1}$, you get the same system as above.
For numerical methods, the second system has the disadvantage, that you can't directly apply most ODE-Solver to that system. So I don't see any reason to write it in this form.
But we can just try to create an Explicit Euler Method for these kind of matrix systems:
Replace $$u_i' \approx \frac{u_i(t+h) - u(t)}{h} = \frac{u_i^{(t+1)} - u^{t}}{h} $$
yields
\begin{align*}
u_1^{(t+1)} &= u_1^t + hF_1 \\
u_2^{(t+1)} &= u_2^t + hF_2 \\
x^2u_1^{(t+1)}+xu_2^{(t+1)}+xu_3^{(t+1)} &= x^2u_1^{(t)}+xu_2^{(t)}+xu_3^{(t)} + hF_3
\end{align*}
And that translates to:
$$Mu^{(t+1)}=Mu^{(t)}+h(\tilde{A}u^{(t)}+\tilde{b}).$$
Again, if you would mutlipy by $M^{-1}$ you would get the Explicit Euler Method for the system $u'=Au+b$.
The reason why this works, is that we assumed $F$ to be a linear function, and the difference quotient is also a linear operator in $f$.
B. In first order, you get the root $T_h$ of a function $f_h(t)=a(t)+hb(t)$ as approximation of a root $T_*$ of $a$. Expanding $f_h(T_*+hv)=a(T_*)+\dot a(T_*)hv+hb(T_*)+O(h^2)$ for $T_h=T_*+hv$ gives an estimate of $v=-\frac{ b(T_*)}{\dot a(T_*)}$.
the numerical approximation $f_h(t)=θ_h(t)$ of step size $h$ seen as $a(t)+hb(t)+...$ with $a(t)=θ_0(t)$ the exact solution, for several values of $h$. It is visible that not only the vertical value has a perturbation proportional to $h$, but also the root location.
So if you know an approximation for $\dot a(T_h)$ and $b(T_h)$, you get $T_h+h\frac{b(T_h)}{\dot a(T_h)}$ as improved root estimation of the root $T_*$ of $a=θ_0$. The important part is that $h\frac{b(T_h)}{\dot a(T_h)}$ is an error estimate of $T_h$.
$\dot a(T_h)=\dot θ_0(T_h)$ you get directly from the differential equation, $b(T_h)$ can be estimated by comparing the results for two different step sizes.
Which raises the question of if it is simpler to just estimating the error of $T_h$ by comparing it to $T_{2h}$. So compute the error for some relatively large but still reasonable $h$ and then scale $h$ so that the expected scaled-down error is in the desired region.
Numerical values for $T_h$ with a secant line. The slope for small $h$ is a little smaller than $0.5$, but still this rough estimate is sufficient to determine $h=10^{-3}$ as sufficient to get 3 correct digits after the dot.
\begin{array}{c|c}
h&T_h\\\hline
0.000500 & 6.70013638\\
0.001000 & 6.70029805\\
0.002000 & 6.70062424
\end{array}
C. just asks for bounds on $|b(T)|$ based on the global error formula $$e(T)\le\frac{M_2}{2L}(e^{LT}-1)h$$ where $M_2$ is a bound for the second derivative around the solution and $L$ the Lipschitz constant. $\dot a(T)$ can again be used directly, if you are looking for a root of $a(t)=\dot θ(t)$, then the value of $\dot a(t)=\ddot θ(t)=-\sinθ(t)$ is known approximately because $θ(T)$ will still be close to the maximum amplitude.
Best Answer
Make a little table -- I've filled in the first couple of rows for you:
\begin{array}{|c|c|c|c|c|} \hline x & y & \Delta x & \frac{dy}{dx} = x - y & \Delta y \approx \frac{dy}{dx}\Delta x \\ \hline 0 & 1 & 0.2 & -1 & -0.2 \\ \hline 0.2 & 0.8 & 0.2 & -0.6 &-0.12 \\ \hline 0.4 & 0.68 & 0.2 & & \\ \hline 0.6 & &0.2 & & \\ \hline 0.8 & &0.2 & & \\ \hline 1 & & 0.2 & & \\ \hline \end{array}
You are done when you get to the bottom left.