For subquestion (1), this site already has How-tos.
For subquestions (2) and (3), here is a scheme.
Suppose you have $n$ mass points that can move along the $x$ axis.
Let $m_1,\ldots,m_n$ denote their (constant) masses.
Let $x_1,\ldots,x_n$ denote the displacements of the mass points
from their rest position.
No springs, no dampers yet.
At this stage, your system of equations looks like this:
$$\underbrace{\begin{pmatrix}m_1\\&m_2\\&&\ddots\\&&&m_n\end{pmatrix}}_{\mathbf{M}}
\underbrace{\begin{pmatrix}\ddot{x}_1\\\ddot{x}_2\\\vdots\\\ddot{x}_n\end{pmatrix}}_{\mathbf{\ddot{x}}}
+\underbrace{\begin{pmatrix}\cdot&\cdot&\cdots&\cdot\\\cdot&\cdot&\cdots&\cdot\\\vdots&\vdots&\ddots&\vdots\\\cdot&\cdot&\cdots&\cdot\end{pmatrix}}_{\mathbf{B}}
\underbrace{\begin{pmatrix}\dot{x}_1\\\dot{x}_2\\\vdots\\\dot{x}_n\end{pmatrix}}_{\mathbf{\dot{x}}}
+\underbrace{\begin{pmatrix}\cdot&\cdot&\cdots&\cdot\\\cdot&\cdot&\cdots&\cdot\\\vdots&\vdots&\ddots&\vdots\\\cdot&\cdot&\cdots&\cdot\end{pmatrix}}_{\mathbf{K}}
\underbrace{\begin{pmatrix}x_1\\x_2\\\vdots\\x_n\end{pmatrix}}_{\mathbf{x}}
= \underbrace{\begin{pmatrix}f_1\\f_2\\\vdots\\f_n\end{pmatrix}}_{\mathbf{f}}
$$
Here, empty matrix entries are to be interpreted as zeros.
$\mathbf{M}$ is the mass matrix (and it's diagonal here),
$\mathbf{B}$ is the damping matrix (empty yet),
$\mathbf{K}$ is the stiffness matrix (empty yet),
$\mathbf{x}$ is the time-dependent vector of the displacements
that you want to describe, and $\mathbf{f}$ is
a time-dependent vector of forces: $f_1$ acts directly on the
first mass point (mass $m_1$), $f_2$ acts directly on the second,
and so on. You may assume $\mathbf{f}$ to be zero; but we will sometimes
think of adding some forces, that's why I have given names to them.
For simplicity, let us suppose that we have only three mass points, $n=3$.
This means that I do not need to write those $\ddots$ dots again.
Now let us add (linear) parts.
Say, you have a spring with stiffness $k_1$,
and that spring shall be mounted between mass point $m_3$ and the wall
(or the ground, whatever, something fixed in any case).
To represent this, we could add to $f_3$ the term $-k_1 x_3$,
but everything that has to do with some $x_i$
shall be moved over to the left-hand side, so we write:
$$\underbrace{\begin{pmatrix}m_1\\&m_2\\&&m_3\end{pmatrix}}_{\mathbf{M}}
\underbrace{\begin{pmatrix}\ddot{x}_1\\\ddot{x}_2\\\ddot{x}_3\end{pmatrix}}_{\mathbf{\ddot{x}}}
+\underbrace{\begin{pmatrix}\cdot&\cdot&\cdot\\\cdot&\cdot&\cdot\\\cdot&\cdot&\cdot\end{pmatrix}}_{\mathbf{B}}
\underbrace{\begin{pmatrix}\dot{x}_1\\\dot{x}_2\\\dot{x}_3\end{pmatrix}}_{\mathbf{\dot{x}}}
+\underbrace{\begin{pmatrix}\cdot&\cdot&\cdot\\\cdot&\cdot&\cdot\\\cdot&\cdot&k_1\end{pmatrix}}_{\mathbf{K}}
\underbrace{\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}}_{\mathbf{x}}
= \underbrace{\begin{pmatrix}f_1\\f_2\\f_3\end{pmatrix}}_{\mathbf{f}}
$$
Note that $k_1$ lands on the main diagonal of $\mathbf{K}$ in the row
(and column) corresponding to $m_3$.
Let us add another part. A linear damper with damping constant $b_1$ shall be mounted between mass points $m_1$ and $m_2$. Again, we could think of this as
adding to $f_1$ the expression $-b_1(\dot{x}_1-\dot{x}_2)$ and adding to $f_2$
the expression $-b_1(\dot{x}_2-\dot{x}_1)$, but that stuff is to be moved over
to the left-hand side, therefore we get:
$$\underbrace{\begin{pmatrix}m_1\\&m_2\\&&m_3\end{pmatrix}}_{\mathbf{M}}
\underbrace{\begin{pmatrix}\ddot{x}_1\\\ddot{x}_2\\\ddot{x}_3\end{pmatrix}}_{\mathbf{\ddot{x}}}
+\underbrace{\begin{pmatrix}b_1&-b_1&\cdot\\-b_1&b_1&\cdot\\\cdot&\cdot&\cdot\end{pmatrix}}_{\mathbf{B}}
\underbrace{\begin{pmatrix}\dot{x}_1\\\dot{x}_2\\\dot{x}_3\end{pmatrix}}_{\mathbf{\dot{x}}}
+\underbrace{\begin{pmatrix}\cdot&\cdot&\cdot\\\cdot&\cdot&\cdot\\\cdot&\cdot&k_1\end{pmatrix}}_{\mathbf{K}}
\underbrace{\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}}_{\mathbf{x}}
= \underbrace{\begin{pmatrix}f_1\\f_2\\f_3\end{pmatrix}}_{\mathbf{f}}
$$
Notice the pattern $\begin{pmatrix}b&-b\\-b&b\end{pmatrix}$ added to the rows and columns of $\mathbf{B}$ that are associated with the mass points
that we connect with the damper. You will see this again.
Let us add more: A linear damper with damping constant $b_2$
between mass point $m_1$ and ground,
and a spring with stiffness $k_2$ between mass points $m_2$ and $m_3$:
$$\underbrace{\begin{pmatrix}m_1\\&m_2\\&&m_3\end{pmatrix}}_{\mathbf{M}}
\underbrace{\begin{pmatrix}\ddot{x}_1\\\ddot{x}_2\\\ddot{x}_3\end{pmatrix}}_{\mathbf{\ddot{x}}}
+\underbrace{\begin{pmatrix}b_1+b_2&-b_1&\cdot\\-b_1&b_1&\cdot\\\cdot&\cdot&\cdot\end{pmatrix}}_{\mathbf{B}}
\underbrace{\begin{pmatrix}\dot{x}_1\\\dot{x}_2\\\dot{x}_3\end{pmatrix}}_{\mathbf{\dot{x}}}
+\underbrace{\begin{pmatrix}\cdot&\cdot&\cdot\\\cdot&k_2&-k_2\\\cdot&-k_2&k_1+k_2\end{pmatrix}}_{\mathbf{K}}
\underbrace{\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}}_{\mathbf{x}}
= \underbrace{\begin{pmatrix}f_1\\f_2\\f_3\end{pmatrix}}_{\mathbf{f}}
$$
Note that the matrix entries are added up when the patterns overlap.
Let us be a bit sophisticated now and add a spring with stiffness $k_3$
between mass points $m_1$ and $m_3$, bypassing $m_2$ somehow. Oh, and
$m_2$ should get a damper $b_3$ with the other side fixed. No problem:
$$\underbrace{\begin{pmatrix}m_1\\&m_2\\&&m_3\end{pmatrix}}_{\mathbf{M}}
\underbrace{\begin{pmatrix}\ddot{x}_1\\\ddot{x}_2\\\ddot{x}_3\end{pmatrix}}_{\mathbf{\ddot{x}}}
+\underbrace{\begin{pmatrix}b_1+b_2&-b_1&\cdot\\-b_1&b_1+b_3&\cdot\\\cdot&\cdot&\cdot\end{pmatrix}}_{\mathbf{B}}
\underbrace{\begin{pmatrix}\dot{x}_1\\\dot{x}_2\\\dot{x}_3\end{pmatrix}}_{\mathbf{\dot{x}}}
+\underbrace{\begin{pmatrix}k_3&\cdot&-k_3\\\cdot&k_2&-k_2\\-k_3&-k_2&k_1+k_2+k_3\end{pmatrix}}_{\mathbf{K}}
\underbrace{\begin{pmatrix}x_1\\x_2\\x_3\end{pmatrix}}_{\mathbf{x}}
= \underbrace{\begin{pmatrix}f_1\\f_2\\f_3\end{pmatrix}}_{\mathbf{f}}
$$
So this is the scheme. Things become more complicated when objects can move in
more than one dimension and when parts can be rotated, but in the end, system matrices will be assembled from element matrices.
You have an ODE
$$
\frac{dx}{dt} = f(t,x) = (x+1)t
$$
which separates to
$$
\frac{dx}{x+1} = tdt\\
\log |x+1| + C = \frac{t^2}{2}\\
x = C_1 e^{t^2/2} - 1.
$$
With $x(0) = 0$ that gives $C_1 = 1$ and
$$
x(t) = e^{t^2/2} - 1.
$$
Computing the first step we get
$$\begin{aligned}
k_1 &= hf(0, 0) = 0\\
k_2 &= hf\left(\frac{h}{2}, \frac{k_1}{2}\right) = h\frac{h}{2}\\
k_3 &= hf\left(\frac{h}{2}, \frac{k_2}{2}\right) = h\frac{h(h^2 + 4)}{8}\\
k_4 &= hf\left(h, k_3\right) =
h\frac{h(8+4h^2+h^4)}{8}\\
x_{1} &= \frac{1}{6}(k_1 + 2k_2 + 2k_3 + k_4) =
\frac{1}{48} h^2 \left(h^4+6
h^2+24\right) = \\
&= \frac{h^2}{2} + \frac{h^4}{8} + O(h^6).
\end{aligned}
$$
And the exact solution is
$$
x(h) = e^{h^2/2} - 1 = \frac{h^2}{2} + \frac{h^4}{8} + O(h^6).
$$
This agrees with the fact that method is of fourth order, thus the local truncation error is $O(h^5)$.
Best Answer
I wanted to note that the method you are using is not foward Euler. What you are doing is:
$$v(t+\delta t) = v(t) + \delta t \cdot a(t) $$ $$x(t+\delta t) = x(t) + \delta t \cdot v(t+\delta t)$$
which is called symplectic Euler method, and has much better convergence properties than the standard Forward Euler method. Forward Euler would be:
$$v(t+\delta t) = v(t) + \delta t \cdot a(t) $$ $$x(t+\delta t) = x(t) + \delta t \cdot v(t)$$
For your answer, the trick is indeed to consider the velocity as a final variable of your system and write what @lhf said:
$$u(t) = (x(t), v(t))$$ $$u' = (v, F/m)$$
You have to solve a single system with six equations per points now. But each equation is only a first order ODE.