Given the vector-valued differential equation \mathbf{y}'(t) = \mathbf{f}(\mathbf{y}(t)), the Euler backward scheme is given by discretizing the time derivative and evaluating $\mathbf{f}$ at the value of $\mathbf{y}$ to be determined:
$\frac{\mathbf{y}_{n+1} - \mathbf{y}_n}{h} = \mathbf{f}(\mathbf{y}_{n+1})$
where $h$ is the time step and $t = nh$, so you have
$\mathbf{y}_{n+1} - \mathbf{y}_n - h\mathbf{f}(\mathbf{y}_{n+1}) = 0$
which is an implicit equation for $\mathbf{y}_{n+1}$. There are two interesting facts about the case you're trying to solve:
(1) the equation is basically linear, i.e. $\mathbf{f}(\mathbf{y}) = \mathbf{b} + A\mathbf{y}$ where $\mathbf{b}=(0,0,1,0)$ and $A$ is a matrix.
(2) the differential equation for $y_3$ and $y_4$ has no dependence on $y_1$ and $y_2$, implying that you can solve for $y_3$ and $y_4$ separately, and then calculate $y_1$ and $y_2$.
If we consider that $\mathbf{f}$ is linear, then we get the equation
$\mathbf{y}_{n+1} - \mathbf{y}_n - hA\mathbf{y}_{n+1} - h\mathbf{b} = 0$
which can be solved explicitly for $\mathbf{y}_{n+1}$, giving
$\mathbf{y}_{n+1} = (I - hA)^{-1}(\mathbf{y}_n + h\mathbf{b})$
where $I$ is the identity matrix. This requires finding the inverse of a 4 x 4 matrix which, whilst it is entirely possible using a computer, isn't necessarily desirable for stability reasons. It would be better to just consider the last two components of the equation, the $y_3$ and $y_4$ components. After manually inverting the 2 x 2 matrix you get
$y_3^{n+1} = \frac{y_3^n + h y_4^n + h}{1+h^2}$
$y_4^{n+1} = \frac{y_4^n - hy_3^n}{1+h^2}$
which you should easily be able to program in Matlab, and then use the values for $y_3$ and $y_4$ to calculate $y_2$ and $y_1$.
Note that I've given you a very specific implementation of backward Euler which is tailored for the problem at hand. A more general backward Euler solver would have to use a root-finding technique such as Newton-Raphson iteration to calculate $\mathbf{y}_{n+1}$ at each step.