Extending numerical Euler method to higher order differential equations 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$$
 A: 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$. 
A: Indeed you would use the first form
$$
xu_3'+xu_3+x^2u_2+u_1=x
$$
so you have exactly one derivative in the equation.
As you only perform one step, you get
$$
u_1(h)= u_1(0)+hu_2(0)
$$
or 
$$
y(0.1)=y(0)+0.1\cdot y'(0)
$$
which can be evaluated only using the intial values without using the equation in any way.
