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)$.
As suggested in the comments, higher-order ODEs are rewritten as a system of first-order ODEs. To illustrate this, an example is provided below.
Let us consider the ODE $x^{(4)} + x^{(3)} + x'' + x' + x = y(t)$, with $x(t_0) = x_0$, $x'(t_0) = \dot x_0$, $x''(t_0) = \ddot x_0$, $x^{(3)}(t_0) = \dddot x_0$. Using the vector notation $\boldsymbol{x} = (x,x',x'',x^{(3)})^\top = (x_1,x_2,x_3,x_4)^\top$, the ODE is rewritten
$$
\underbrace{
\left(
\begin{array}{c}
x'_1 \\
x'_2 \\
x'_3 \\
x'_4
\end{array}
\right)}_{\boldsymbol{x}'}
=
\underbrace{
\left(
\begin{array}{c}
x_2 \\
x_3 \\
x_4 \\
y(t) - x_4 - x_3 - x_2 - x_1
\end{array}
\right)}_{\boldsymbol{f}(\boldsymbol{x},t)} ,
$$
with the initial condition $\boldsymbol{x}(t_0) = (x_0,\dot x_0,\ddot x_0, \dddot x_0)^\top = \boldsymbol{x}_0$. Then the Runge-Kutta 4 method writes
$$
\boldsymbol{x}_{n+1} = \boldsymbol{x}_{n} + \frac{h}{6} \left(\boldsymbol{k}_1 + 2\boldsymbol{k}_2 + 2\boldsymbol{k}_3 + \boldsymbol{k}_4\right),
$$
where
$$
\begin{aligned}
\boldsymbol{k}_1 &= \boldsymbol{f}(\boldsymbol{x}_n, t_n) \, ,\\
\boldsymbol{k}_2 &= \boldsymbol{f}(\boldsymbol{x}_n + \tfrac{h}{2}\boldsymbol{k}_1, t_n + \tfrac{h}{2}) \, ,\\
\boldsymbol{k}_3 &= \boldsymbol{f}(\boldsymbol{x}_n + \tfrac{h}{2}\boldsymbol{k}_2, t_n + \tfrac{h}{2}) \, ,\\
\boldsymbol{k}_4 &= \boldsymbol{f}(\boldsymbol{x}_n + h \boldsymbol{k}_3, t_n + h ) \, .
\end{aligned}
$$
Best Answer
You have $f(\binom{x}{y})=\binom{y^2}{x^2+xy}$. Give the state vector a different name, $v=\binom{x}{y}$. Then $$ k_1=hf(v_0) $$
The vector $k_2$ is computed at the position $v_0+\frac12k_1$, $$ k_2=hf(v_0+\tfrac12k_1)=\pmatrix{h(y_0+\tfrac12k_{y1})^2\\h(x_0+\tfrac12k_{x1})^2+h(x_0+\tfrac12k_{x1})(y_0+\tfrac12k_{y1})} $$ etc.
It is useful to have the evaluation of $f$ as a separate function so that one is forced to compute the intermediate positions only once.