How to solve a 4th order differential equation in space and time variables in Julia

I am relatively new to solving differential equations in Julia and thus cant figure out how to solve a higher order ode in 2 independent variables(space and time) and thus would request for assistance.

I am trying to plot the curve using Julia between h(height of fluid film) vs x(length span) of a fluid film suspended from a ceiling that drains into drops due to the Rayleigh–Taylor instability. There are 2 key differential equations that govern this phenomena. The two odes are listed here

Here h represents the height of the fluid film ,B0=0.134 , and the rest are derivatives of h and q with respect to time(t) and space(x).The notation Dxxxh means 3rd order derivative of h wrt to x. The space span can be considered as L=24. The boundary conditions are defined here.

The initial value of h at t=0 can be found using this expression where epsilon=0.0009.The expected plot is as shown here. The plot slides as time progresses.

NeuralPDE.jl is a physics-informed neural network based PDE solver. It can handle this kind of equation when given symbolically. Here is an example solving the Kuramoto–Sivashinsky equation, also a 4th order PDE.

using NeuralPDE, Flux, ModelingToolkit, GalacticOptim, Optim, DiffEqFlux
import ModelingToolkit: Interval, infimum, supremum

@parameters x, t
@variables u(..)
Dt = Differential(t)
Dx = Differential(x)
Dx2 = Differential(x)^2
Dx3 = Differential(x)^3
Dx4 = Differential(x)^4

α = 1
β = 4
γ = 1
eq = Dt(u(x,t)) + u(x,t)*Dx(u(x,t)) + α*Dx2(u(x,t)) + β*Dx3(u(x,t)) + γ*Dx4(u(x,t)) ~ 0

u_analytic(x,t;z = -x/2+t) = 11 + 15*tanh(z) -15*tanh(z)^2 - 15*tanh(z)^3
du(x,t;z = -x/2+t) = 15/2*(tanh(z) + 1)*(3*tanh(z) - 1)*sech(z)^2

bcs = [u(x,0) ~ u_analytic(x,0),
       u(-10,t) ~ u_analytic(-10,t),
       u(10,t) ~ u_analytic(10,t),
       Dx(u(-10,t)) ~ du(-10,t),
       Dx(u(10,t)) ~ du(10,t)]

# Space and time domains
domains = [x ∈ Interval(-10.0,10.0),
           t ∈ Interval(0.0,1.0)]
# Discretization
dx = 0.4; dt = 0.2

# Neural network
chain = FastChain(FastDense(2,12,Flux.σ),FastDense(12,12,Flux.σ),FastDense(12,1))

discretization = PhysicsInformedNN(chain, GridTraining([dx,dt]))
@named pde_system = PDESystem(eq,bcs,domains,[x,t],[u(x, t)])
prob = discretize(pde_system,discretization)

cb = function (p,l)
    println("Current loss is: $l")
    return false
end

opt = Optim.BFGS()
res = GalacticOptim.solve(prob,opt; cb = cb, maxiters=2000)
phi = discretization.phi