There is absolutely nothing wrong with converting the second-order system
to first-order form and then using appropriate numerical methods to solve it.
Both implicit and explicit Euler methods can be used. However, both have only
first-order accuracy, i.e. the error is proportional to the time step size. And, of
course, explicit Euler is not stable unless the time step is sufficiently small.
A method with second-order accuracy can be obtained by writing the solution of
$$
\frac{dy}{dt} = f(y,t)
$$
as
$$
y_{n+1} = y_n + h [\theta ( t_n, y_n) + (1 - \theta ) f ( t_{n+1}, y_{n+1})]
$$
where $h$ is the time step, $h=t_{n+1}-t_n$.
Selecting $\theta=0$ gives implicit Euler, $\theta=1$ gives explicit Euler, and
$\theta=1/2$ gives the second-order accurate trapezoidal method.
The Newmark numerical "method" for solving second-order ODE is actually a family
of methods where a specific method can be chosen by selecting values for the $\beta$
and $\gamma$ parameters. One of the more popular methods in this family is obtained by
selecting $\beta=1/4$ and $\gamma=1/2$ which gives a constant acceleration and a linear
velocity over the time step. Using this method on the second-order equations gives identical results to using the
trapezoidal method on the first order version of the equations.
I would not say the that elasticity equations are inherently stiff. The stiffness results
when they are discretized in space using FEM or some other approach. The system of
ODE becomes stiffer as the spatial mesh is refined. This is why implicit ODE methods
are often used for these types of problems.
