If it was a second order DE I can follow the logic however I get very confused when trying to numerically solve a third order DE. I have to estimate $y(0.1)$ using step size 0.1
$$xy''' + xy'' + x^2y' + y = x,\quad with \quad y(0) = 1, y'(0) = 0, y''(0) = 1$$
where
$$y_0 = y(x_0)$$ $$y_{n +1} = y_n + hf(x_n, y_n)$$
$$x_{n +1} = x_n + h$$
Now I can set $u_1=y, u_2=y', u_3=y'' $ which gives $u_1'=y'=u_2$ and $u_2'=y''=u_3$ and $u_3'=y'''$
So I can sub these in to get
$$xu_3'+xu_3+x^2u_2+y=x\tag{*}$$ but I can also get
$$xu_3'+xu_2'+x^2u_1'+y=x\tag{**}$$
Would they both lead to the same answers? The first one does seem easier to work with.
And assuming I go with the first one, I still do not know how to set up my iteration because I have no idea how to relate them.
Approximate the value of $y (0.2)$ , where $y(x)$ is the solution of the
initial-value problem
$$y''+ xy' + y = 0, y(0) = 1, y'(0) = 2$$
Substitute $u = y'$ to obtain the system
$y' = u, u ' = −xu − y$
Euler’s scheme then gives
$$y_{n +1}= y_n + hu_n$$
$$u_{n +1}= u_n + h( −x_nu_n − y_n)$$
If we use the step-size $h = 0.1$ and $y_0 = 1$, $u_0= 2$, we obtain
$$y_1 = y_0 + (0.1)u_0 = 1 + (0.1)2 = 1.2$$
$$u_1 = u_0 + (0.1)( −x 0u 0 − y_0) = 2 + (0.1)( −1) = 1.9$$
$$y_2 = y_1 + (0.1)u_1 = 1.2 + (0.1)(1.9) = 1.39$$
$$u_2 = u_1 + (0.1)( −x_1 u_1 − y_1) = 1.9 − (0.1)1.39 = 1.761.1.9 + (0.1)$$
Thus we obtain $$y(0.2) ≈ 1.39$$ and $$y'(0.2) ≈ 1.761$$
Best Answer
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$.