Error in program for coupled PDE

Hello
Could somebody tell me what went wrong? Here the program

using OrdinaryDiffEq, ModelingToolkit, MethodOfLines, DomainSets, Plots,
Optimization, OptimizationPolyalgorithms, SciMLSensitivity
ForwardDiff,JLD2

# Parameters, variables, and derivatives
@parameters t x
@variables u(..)  v(..) ud(..) vd(..)
xmin = 0.1
xmax = 1.0
Dt = Differential(t)
Dxx = Differential(x)^2
Dx = Differential(x)
# 1D PDEs and boundary conditions
eq  = [Dt(u(t, x)) ~ ud(t, x), Dt(v(t, x)) ~ vd(t, x),
       Dt(ud(t, x)) ~ v(t, x)^4*Dxx(u(t, x))-5*v(t, x)^4*Dx(u(t, x))^2/(3*u(t, x))
       +5*ud(t, x)^2/(3*u(t, x)),
       Dt(vd(t, x)) ~ v(t, x)^4*(Dxx(v(t, x))+Dx(v(t, x))^2/v(t, x)+3*Dx(v(t, x)))+3*vd(t, x)^2/v(t, x)
+(2/u(t, x))*(v(t, x)^5*(Dxx(u(t, x)))-(5/3)*ud(t, x)^2/(u(t, x))+Dx(u(t, x))                                                        
+v(t, x)^4*Dx(u(t, x))*Dx(v(t, x))-ud(t, x)*vd(t, x)
)]                                                                               
bcs = [u(0, x) ~ 1/(x+2),
       v(0, x) ~ x^5,
       ud(0, x) ~ 0,
       vd(0, x) ~ 0, 
        u(t, xmin) ~ 1/((xmin+2)*(t+2)),
        v(t, xmin) ~ xmin^5/(t+2)^4,
        u(t, xmax) ~ 1/((xmax+2)*(t+2)),
        v(t, xmax) ~ xmax^5/(t+2)^4,
        ud(t, xmin) ~ 0,
        vd(t, xmin) ~ 0,
        ud(t, xmax) ~ 0,
        vd(t, xmax) ~ 0]

# Space and time domains
domains = [t ∈ Interval(0.0, 1.0),x ∈ Interval(xmin, xmax)]
# PDE system
@named pdesys = PDESystem(eq, bcs, domains, [t, x], [u(t, x),v(t, x),ud(t, x),vd(t, x)])
# Method of lines discretization
dx = 0.1
order = 4
discretization = MOLFiniteDifference([x => dx], t)
# Convert the PDE problem into an ODE problem
prob = discretize(pdesys,discretization)
# Solve ODE problem
using OrdinaryDiffEq
sol = solve(prob, Tsit5(), saveat=0.2)
# Plot results 
discrete_x = sol[x]
discrete_t = sol[t]
solu = sol[u(t, x)]
solv = sol[v(t, x)]
x = range(1, length=11)
t = range(1, length=6)
p1=plot(x, t, solu, st=:surface, title="3D Surface Plot", xlabel="X", ylabel="t", zlabel="u")
p2=plot(x, t, solv, st=:surface, title="3D Surface Plot", xlabel="X", ylabel="t", zlabel="v")
plot(p1,p2)

#savefig("HeatDirich2.png")                                                   

It has to do with allocation of the matrices