solution of system of coupled partial differential equations

Answer 1

[Precise Simulation]

One way to approach this is, due to the 2nd derivative in time, to split the time dependent variables into

du1/dt = u1h
du1h/du = rest of the equation

as typical for the wave equation fem simulations, and using Hermite finite element shape functions to account for the 4th order derivative, as done here for the Euler-Bernoulli beam deformation simulation.