Method Error when using Method of Lines to solve PDEs

Hello all! I’m new to Julia and am trying to use the Method of Lines and MTK packages to solve a system of PDEs. I am running into repeated method errors, and the stacktrace keeps pointing to a line from Macro Expanded code, and I’m not sure where to start with troubleshooting. Here’s the code and stack trace.

using ModelingToolkit
using DifferentialEquations
using MethodOfLines
using DomainSets
using Symbolics

#Parameters and Variables:
JJ = 3
N = 45
P = 20
@parameters r z t
@parameters D ε εp K kI u_avg qm Cf rp[1:JJ] psd[1:JJ] kf[1:JJ] zL
@variables c(..)[1:JJ] q(..)[1:JJ] Csl(..) Csh(..) #Pore Fluid Concentration, Sorption Concentration, Bulk Concentration of Slice, Bulk Concentration of Column
@variables c_surface(..)[1:JJ] #Helper Variable for cj(r=rpj, t)
@variables Q(..) #Total Pore Concentration 
@variables cstar(..)[1:JJ] #Equilibrium Concentration

#Differential Operators:
Dt = Differential(t)
Dr = Differential(r)
Drr = Differential(r)^2
Dz = Differential(z)

#Equations
eqs = [] 
bcs = []
rp_vals = [50,62,75]
kf_vals = 1.54e-5./rp_vals
frac = [0.22, 0.54, 0.24]
psd_vals = kf_vals./rp_vals .*frac

push!(eqs, Dt(εp *c(r,z,t)[j]) + Dt(q(r,z,t)[j]) ~ D*(2/r*Dr(c(r,z,t)[j]) + Drr(c(r,z,t)[j])) for j in 1:JJ ) #Diffusion into a pore
push!(eqs, c_surface(z, t)[j] ~ c(rp[j],z, t)[j] for j in 1:JJ) #Surface Flux
push!(eqs, cstar(r,z,t)[j] ~ q(r,z,t)[j] / (K * (qm - q(r,z,t)[j])) for j in 1:JJ) #Langmuir Isotherm Equilibrium Relationship
push!(eqs, Dt(q(r,z,t)[j] ~ kI * (c_surface(z,t)[j] - cstar(r,z,t)[j])) for j in 1:JJ) #LDF Model for Equilibrium
push!(eqs, Q(z,t) ~ (εp * c(r,z,t)[j] + q(r,z,t)[j]) for j in 1:JJ) #Total Concentration entering Pore 


push!(eqs, Dt(Q(z,t) * -(1/3) * (1-ε / ε)) ~ sum(psd[j] * (Csl(z,t) - c_surface(z,t)[j]) for j in 1:JJ) ) # Column Material Balance
push!(eqs, Dt(Csh(z,t)) ~ (-(1-ε)/ε) * Dt(Q(z,t)) - (u_avg / ε) * Dz(Csl(z,t))) #Column Balance

#Boundary Conditions

push!(bcs, Dr(c(0, z, t)[j]) ~ 0 for j in JJ) #Inside of Particle Neumann 
push!(bcs, Dr(c(rp[j],z, t)[j]) ~ kf[j] * (Csl(z, t) - c_surface(z,t)[j]) for j in JJ ) #Surface of Particle Robin

push!(bcs, Dz(Csl(zL, t)) ~ 0) #Bottom of Column Neumann
push!(bcs, Csl(0, 0) ~ Cf) #Feed Concentration at top of column Dirichlet at both z = 0, t = 0

#Space and Time Domains
dmns = []

push!(dmns, r ∈ Interval(0, maximum(rp_vals))) #Interval for particle radii
push!(dmns, z ∈ Interval(0, zL)) #Interval for axial dimension from 0 to column length zL
push!(dmns, t ∈ Interval(0, Inf)) #Time Interval, based on spatial bounds


#Defining the PDE Systems
@named pdesys = PDESystem(eqs, bcs, dmns, [r,z,t], [c(r,z,t)[1:JJ], q(r,z,t)[1:JJ], Csl(z,t), Csh(z,t), Q(z,t), c_surface(z,t)[1:JJ], cstar(r,z,t)[1:JJ]], 
[D => 5e-6, ε => 0.25, εp => 0.1, K => 367, kI => 4000, u_avg => 2.5, qm => 242, Cf => 3.0, rp[1:JJ] .=> rp_vals, psd[1:JJ] .=> psd_vals, kf[1:JJ] .=> kf_vals, zL => 2 ])

#Method of Lines Discretizations
order = 2
discretization = MOLFiniteDifference([r => N, z => P], t, approx_order = order)

#Convert the Problem from a PDE to an ODE
@time problem = discretize(pdesys, discretization)

#Solve the ODE
using DifferentialEquations
@time solution = solve(problem, FBDF(), saveat = 0.1)

#Plot Final Concentration vs Time
using Plots
plot(solution[Csh(z,t)], linewidth = 2, xaxis = "Time(t)", yaxis = "Concentration")

and here’s the stack trace

ERROR: MethodError: no method matching real(::Base.Generator{UnitRange{Int64}, var"#11#12"})
The function `real` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  real(::Type{Union{}}, Any...)
   @ Base complex.jl:123
  real(::Symbolics.SemiMonomial)
   @ Symbolics C:\Users\U1087908\.julia\packages\Symbolics\lJiB2\src\semipoly.jl:123
  real(::Missing)
   @ Base missing.jl:101
  ...

Stacktrace:
  [1] hascomplex(term::Base.Generator{UnitRange{Int64}, var"#11#12"})
    @ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\symbolic_utils.jl:292
  [2] (::PDEBase.var"#125#133")(eq::Base.Generator{UnitRange{Int64}, var"#11#12"})
    @ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\make_pdesys_compatible.jl:100
  [3] _any(f::PDEBase.var"#125#133", itr::Vector{Any}, ::Colon)
    @ Base .\reduce.jl:1237
  [4] any(f::Function, a::Vector{Any}; dims::Function)
    @ Base .\reducedim.jl:992
  [5] handle_complex(pdesys::PDESystem)
    @ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\make_pdesys_compatible.jl:100
  [6] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
    @ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\symbolic_discretize.jl:10
  [7] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, kwargs::@Kwargs{})
    @ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\discretization_state.jl:57
  [8] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
    @ PDEBase C:\Users\U1087908\.julia\packages\PDEBase\CpWUp\src\discretization_state.jl:54
  [9] macro expansion
    @ .\timing.jl:581 [inlined]
 [10] top-level scope
    @ c:\Users\U1087908\OneDrive - CompanyName\Vscode\ChromatographyModeling\BreakthroughPredict.jl:315
    @ Base .\reduce.jl:1237

my guess is it’s an issue with the syntax for the indexing, but I’m not sure how to troubleshoot that. Thank you for any help.