# Issue with (2+1)D Diffusion Reaction PDE with DiffEqOperators.jl

**URL:** <https://discourse.julialang.org/t/issue-with-2-1-d-diffusion-reaction-pde-with-diffeqoperators-jl/74840>\
**Category:** Numerics\
**Tags:** differentialequation\
**Created:** [January 19, 2022, 1:25am UTC](https://discourse.julialang.org/t/issue-with-2-1-d-diffusion-reaction-pde-with-diffeqoperators-jl/74840 "2022-01-19T01:25:47Z")\
**Posts on this page:** 1\
**Page:** 1

<div class="post-metadata">

**Author:** ![ash](https://avatars.discourse-cdn.com/v4/letter/a/f19dbf/32.png) [@ash](https://discourse.julialang.org/u/ash)\
**Post date:** [January 19, 2022, 1:25am UTC](https://discourse.julialang.org/t/issue-with-2-1-d-diffusion-reaction-pde-with-diffeqoperators-jl/74840/1 "2022-01-19T01:25:47Z")

</div>

I am studying a (2+1)D diffusion-reaction PDE. I have asked about this before and got plenty of help and that resolved my issues. However, I just realized I was getting the wrong results because the model should be in two spatial dimensions, not one, so I started working on modifying it.

Here’s a minimal non-working example, with a much simpler PDE:

```julia
using DifferentialEquations, DiffEqOperators, LinearAlgebra, Symbolics, SparseArrays

# Define FD operator
N = 60
δ = 30.0/N
x = range(0, step=δ, length=N)
y = x
ord_deriv = 2
ord_approx = 4
Δx = CenteredDifference{1}(ord_deriv, ord_approx, δ, N)
Δy = CenteredDifference{2}(ord_deriv, ord_approx, δ, N)
bc = Neumann0BC(δ, 2)
Qx, Qy = MultiDimBC(bc, size(rand(N,N)))

# Define ODE
function myModel!(du,u,p,t)
    c₁ = u[1,:,:]
    c₂ = u[2,:,:]
    du[1,:,:] = - c₂.^2 .+ 2.0*(Δx*Qx*c₁ + Δy*Qy*c₁)
    du[2,:,:] = - c₁ .+ 3.0*(Δx*Qx*c₂ + Δy*Qy*c₂)
end

# Prepare the ODE integrator
tspan = (0.0,0.03)
u0 = rand(Float64, (2, N, N))

# Prepare the Jacobian stuff
du0 = copy(u0)
jac_sparsity = Symbolics.jacobian_sparsity((du,u)->myModel!(du,u,p,0.0),du0,u0)
f = ODEFunction(myModel!;jac_prototype=float.(jac_sparsity))
prob = ODEProblem(f,u0,tspan,p)

# Solve!
sol = solve(prob,TRBDF2(linsolve = KLUFactorization()), reltol=1e-8, abstol=1e-10)

```

which leads to the following error:

```julia
MethodError: no method matching *(::MultiDimDirectionalBC{Float64, RobinBC{Float64, Vector{Float64}}, 1, 2, 1}, ::Matrix{Num})

```

which comes from 3rd line in `myModel!`. It seems like a simple issue, it probably has something to do with how I’m indexing `u` but I can’t quite figure it out. Any advice would be appreciated.

If I just do:

```julia
prob = ODEProblem(myModel!,u0,tspan,p)

```

and solve it (i.e., without the fancy Jacobian stuff), it just takes forever or crashes (with the actual PDE under study, not this minimal one).