Solving coupled PDEs with implicit time-dependence

I’m trying to implement a time-dependent Landau model for ferroelectric materials on a 2D cartesian domain, where the time-dependent Landau equation
\frac{dP}{dt} = -\eta (\alpha P+\beta P^3 + \gamma P^5 - k_x \frac{\partial^2 P}{\partial x^2} - k_z \frac{\partial^2 P}{\partial z^2} + P\frac{dV}{dz}) is coupled with the Poisson equation
\nabla^2 V = \frac{1}{\epsilon} \frac{\partial P}{\partial z}
In this system, V is not explicitly dependent on time, and the typical solution routine solves for V for a given P, then uses that solution for V to time step P, and repeats.

I’m wondering which Julia package would be most suited for this task. My first attempt was with MethodOfLines.jl, but the following code froze upon execution (>5min on the discretize step, no output). Any help is appreciated!

mutable struct Ferroelectric
    # Parameters
    ϵ::Float64 # Permittivity
    
    # Landau coefficients
    α::Float64
    β::Float64
    γ::Float64

    # Domain wall coupling coefficients
    k_x::Float64
    k_z::Float64

    # dynamic coefficient
    η::Float64
end

using MethodOfLines, ModelingToolkit, DomainSets, Plots, OrdinaryDiffEq

@parameters t x z
@variables P(..) V(..)

Dt= Differential(t); Dz = Differential(z); Dx = Differential(x)
Dxx = Differential(x)^2; Dzz = Differential(z)^2; Dtt = Differential(t)^2


function simulate(fe::Ferroelectric, lims::Tuple{Float64, Float64, Float64}, ΔV::Float64, Δx::Float64, Δz::Float64, Δt::Float64, order::Int)
    x_lim, z_lim, t_lim = lims

    # Define the PDE System
    eqs = [Dt(P(x, z, t)) ~ -fe.η * (fe.α * P(x, z, t) + fe.β * P(x, z, t)^3 + fe.γ * P(x, z, t)^5 - fe.k_x * Dxx(P(x, z, t)) - fe.k_z * Dzz(P(x, z, t)) + P(x, z, t) * Dz(V(x, z, t))),
    Dxx(V(x, z, t)) + Dzz(V(x, z, t)) ~ P(x, z, t) / fe.ϵ]
    print("System defined \n")

    
    # Define the domain
    domain = [x ∈ Interval(0.0, x_lim), 
              z ∈ Interval(0.0, z_lim),
              t ∈ Interval(0.0, t_lim)]

    # boundary conditions
    bcs = [P(0, z, t) ~ P(x_lim, z, t), # PBC for P
           P(x, z, 0) ~ 5, # initial polarization
           V(0, z, t) ~ V(x_lim, z, t), # PBC for V
           V(x, 0, t) ~ 0,
           V(x, z_lim, t) ~ ΔV, # voltage applied at the top
           V(x, z, 0) ~ ΔV * (z / z_lim), # initial potential distribution
           ]

    @named pdesys = PDESystem(eqs, bcs, domain, [x, z, t], [P(x, z, t), V(x, z, t)])

    # Discretize the PDE System
    discretization = MOLFiniteDifference([x => Δx, z => Δz], t, approx_order=order)

    prob = discretize(pdesys, discretization)
    print("System discretized \n")

    # Solve the PDE System
    sol = solve(prob, Tsit5(), saveat=Δt)

    return sol
end

test_fe = Ferroelectric(16., 1., 1., 1., 0.1, 0.1, 0.1) #arbitrary test values
simulate(test_fe, (1.0, 1.0, 1.0), 1.0, 0.1, 0.1, 0.1, 2)