Write your equations in the form M*DV = RHS where:
and M is the mass matrix. Since your second equation expands to:
the second row of your mass matrix will be [y, x, 0, 0]. Multiply that vector times DV and you'll see that you've recreated the left side of the second differential equation. Generate the remaining rows of the mass matrix similarly.
Create an options structure that specifies the mass matrix using odeset and the 'Mass' name-value pair. The value for that pair will be a function handle that accepts t (time) and the vector V and returns the mass matrix M.
Then write the function that evaluates RHS as a function of t and V. Call the ODE solver specifying that function and the options structure (so the solver knows how to create the mass matrix.