I am trying to solve a system of coupled PDEs shown below:

and

.

I have written some simple code using DifferentialEquations.jl and DiffEqOperators.jl hoping to solve this system:

```
using DifferentialEquations
using DiffEqOperators
# 300x300 grid, x \in 0,30 and t \in 0,0.3
# Define FD operator
Nₓ = 300
Δx = 30.0/(Nₓ)
x = range(0, step=Δx, length=Nₓ)
ord_deriv = 2
ord_approx = 2
Δ = CenteredDifference(ord_deriv, ord_approx, Δx, Nₓ)
bc = Dirichlet0BC(Float64)
#Parameters
f = 1.2
m = 190
q = 0.001
ϵ₁ = 0.01
ϵ = 0.01
K = 1000
Θ = 0.0025
d₀ = 0.01
D₁ = 1
p = [f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ]
function myModel!(du,u,p,t)
f, m, q, ϵ, ϵ₁, D₁, d₀, K, Θ = p
c₁, c₂ = u
du[1] = 1/ϵ*(f*c₂.*(q .- c₁)./(q .+ c₁) .+ c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₁.^2) .+
D₁*Δ*bc*c₁
du[2] = c₁ .* (1 .- m*c₂)./(1 .- m*c₂ .+ ϵ₁) .- c₂ .+ D₁/d₀*(1 + (K-1)/(1+3/Θ))^(-3/2).*Δ*bc*c₂
end
# Prepare the ODE intnegrator
tspan = (0.0,0.3)
u0 = rand(Float64, (1, Nₓ))
prob = ODEProblem(myModel!, u0, tspan, p)
dt = 0.001
solve(prob, Tsit5(), dt = dt)
```

However, I am greeted with this error:

```
ERROR: LoadError: MethodError: no method matching iterate(::DerivativeOperator{Float64, 1, false, Float64, StaticArrays.SVector{3, Float64}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, StaticArrays.SVector{0, StaticArrays.SVector{4, Float64}}, Vector{Float64}, Int64})
```

which makes me think I’m doing something wrong with DiffEqOperators, but it seems to work when you just have one equation. Any thoughts?