# ModelingToolkit - MethodError on isequal()

My apologies if this is something really simple. I have looked at the docs and whatever examples I could find, but I am beginning to confuse the trees for the forest. Any directions would be greatly appreciated.

I am trying to learn ModelingToolkit by way of building a simple, isothermal axial dispersion reactor model with one reaction and three species. I started off by trying a dynamic model, but it looks like the PDE discretization tools are still WIP, so I changed this to a steady-state model.

``````using ModelingToolkit

#=
Differential molar balance for i
∂Cᵢ/∂t = 0 = Dₐ*∂²Cᵢ/∂z² - ∂(uCᵢ)/∂z + ∑(reaction j)νᵢⱼrᵢ
Peₐ = ul/Dₐ (l = characteristic length)

Dankwerts boundary conditions
z = 0   uₛCᵢ(z = 0) - Dₐ*dCₜₒₜ/dz(z=0) = uₛCᵢ(feed)
z = L   dCᵢ/dz(z=L) = 0

Cₜₒₜ = P/RT
Cᵢ = nᵢ/nₜₒₜ*Cₜₒₜ

Reaction:
A + 2B -> C
r = k*pA*pC / (1 + Ka*pA + Kb*pB)
rA = -r
rB = -2r
rC = +r
=#

const N = 3 # Number of species
const ν = [-1.0, -2.0, 1.0]
const R₀ = 8.314
const Cfeed = [1.0, 1.0, 0.0]
const L = 2.0

@parameters T Dₐ k Ka Kb
@variables z u(z) C[1:N](z) p[1:N](z) r[1:N](z)

Dz = Differential(z)
Dzz = Differential(z)^2

# 1D ADPFR with Dankwerts boundary conditions

# Cᵢ = Pᵢ/R₀T
presconc = [C[i] ~ p[i]/R₀/T for i in 1:N]

# A + 2B -> C
# r = k*pA*pC / (1 + Ka*pA + Kb*pB)
reactions = [r[i] ~ ν[i] * (k * p * p / (1 + Ka*p + Kb*p)) for i in 1:N]

# ∂Cᵢ/∂t = 0 = Dₐ*∂²Cᵢ/∂z² - ∂(uCᵢ)/∂z + ∑(reaction j)νᵢⱼrᵢ
axdisp = [0 ~ Dₐ*Dzz(C[i]) - Dz(u*C[i]) + r[i] for i in 1:N]

eqs = vcat(presconc, reactions, axdisp)

sys = ODESystem(eqs)
``````

Each line executes as expected, up to the last, which returns an error:

``````ERROR: LoadError: MethodError: isequal(::Missing, ::Sym{Real, Base.ImmutableDict{DataType, Any}}) is ambiguous. Candidates:
isequal(x, ::SymbolicUtils.Symbolic) in SymbolicUtils at C:\Users\Braam\.julia\packages\SymbolicUtils\4n4Gz\src\types.jl:123
isequal(::Missing, ::Any) in Base at missing.jl:81
Possible fix, define
isequal(::Missing, ::SymbolicUtils.Symbolic)
Stacktrace:
 collect_var!(states::OrderedCollections.OrderedSet{Any}, parameters::OrderedCollections.OrderedSet{Any}, var::Missing, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
@ ModelingToolkit C:\Users\Braam\.julia\packages\ModelingToolkit\3q7jk\src\utils.jl:303
 collect_vars!(states::OrderedCollections.OrderedSet{Any}, parameters::OrderedCollections.OrderedSet{Any}, expr::SymbolicUtils.Add{Real, Int64, Dict{Any, Number}, Nothing}, iv::Sym{Real, Base.ImmutableDict{DataType, Any}})
@ ModelingToolkit C:\Users\Braam\.julia\packages\ModelingToolkit\3q7jk\src\utils.jl:296
 ODESystem(eqs::Vector{Equation}, iv::Nothing; kwargs::Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}})
@ ModelingToolkit C:\Users\Braam\.julia\packages\ModelingToolkit\3q7jk\src\systems\diffeqs\odesystem.jl:172
 ODESystem(eqs::Vector{Equation}, iv::Nothing) (repeats 2 times)
@ ModelingToolkit C:\Users\Braam\.julia\packages\ModelingToolkit\3q7jk\src\systems\diffeqs\odesystem.jl:151
 top-level scope
I think the problem is the term `Dz(u*C[i])`, `ODESystem` works better with simple differential terms. This seems to prevent the error:
``````expanded_eqs = expand_derivatives.(eqs)