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.
