I have played around with your proposed system of equations and solved it with RK45 in scipy. However, the number of function evaluations and the size of the time steps seem be the same for the joint system as for just the $y(t)$ equation. From your explanation I would have expected the joint system to require more evaluations / smaller time steps.

Why can scalar ODEs not be stiff? I thought that $x_t = -\lambda x$ was the prototypical stiff equation where explicit methods needs quite small time steps to remain stable.

Thank you very much for the detailed answer! One thing I did not mention in the question is that I discretized with finite elements in space with triangular elements. Then I get two the mass and stiffness matrices and the equation $Au_{tt} = Bu^{(t)}$ that can be solved for the second t-derivates which I then give to the ODE solver. And here I just used the $A$ and $B$ that came out of the computation. However, my suspicion is that the entries for the boundary points are wrong because they are not surrounded by triangles on all sides. They are on the boundary after all.