# differentialEquation or diffEqOperator, whose problem?

basically, I would like to simulate a gas-solid reaction. u[:,1] is gas reactant, u[:,2] is gas product, and u[:3] is solid product. some times the result is strange.
the code is below

``````using OrdinaryDiffEq, DiffEqOperators, Plots
pyplot()

nknots = 100  # 100 knots
Lbed = 0.03 # bed length
h = Lbed/(nknots+1) #
knots = collect(range(h, step=h, length=nknots))

ord_approx = 2
usg = 0.1
Deg = 2.2E-4

# Δ1, first derivative, (1st order derivative, 2nd order approxi.)
Δ1 = UpwindDifference(1, ord_approx, h, nknots, -usg)
Δ2 = CenteredDifference(2, ord_approx, h, nknots)

bc1 = RobinBC((1.0, -Deg/usg, 1.0), (0.0, 1.0, 0.0), h, ord_approx)
bc2 = RobinBC((1.0, -Deg/usg, 0.0), (0.0, 1.0, 0.0), h, ord_approx)
function flow!(du, u,p,t)
du[:,1] = Δ1*bc1*u[:,1] + Deg*Δ2*bc1*u[:,1]
du[:,2] = Δ1*bc2*u[:,2] + Deg*Δ2*bc2*u[:,2]
end
function source_term!(du,u,p,t)
du .= 0.
for ii = 1:length(du[:,1])
if  u[ii,3] <= 1.0
rA = -0.5*u[ii,1]*(1-u[ii,3])
else
rA = 0.
end
du[ii,1] = rA
du[ii,2] = -rA
du[ii,3] = -rA
end
end
u0 = zeros(length(knots),3)
u0[:,1] .= 0.01
``````

just by running the code below, different results were obtained

``````probSplit = SplitODEProblem(flow!, source_term!, u0, (t0, 4.))
sol_split = solve(probSplit, KenCarp4(autodiff=:false))
plt1 = surface(sol_split.t, knots, Array(sol_split)[:,1,:], xlabel="t [s]", ylabel="L ", camera=(-45, 30), legend=:false)
plt2 = surface(sol_split.t, knots, Array(sol_split)[:,2,:], xlabel="t [s]", ylabel="L ", camera=(-45, 30), legend=:false)
plt3 = surface(sol_split.t, knots, Array(sol_split)[:,3,:], xlabel="t [s]", ylabel="L ", camera=(-45, 30), legend=:false)
plot(plt1, plt2, plt3, layout=(1,3), size = (900, 300))
``````

a normal result

three abnormal results

error

Can anyone tell me where is the problem?

I do not know if it makes sense but `flow!` is missing an assignment for `du[:,3]`.

Good point! They are actually zero.
I thought the `differentialEquations.jl` package will intialize du as zero. But seems not this case. When I add `@. du[:,3] = 0. ` in `flow!` function, everyting is alright.

I think it has changed to make it not work (so the `NaN`s) so that if someone accidentally forgets to assign values, it should not silently work nonetheless.

Have you seen MethodOfLines.jl? It is the successor library to DiffEqOps, and is currently under active development. It already contains most of DEO’s functionality but with much greater flexibility, I recommend that you check it out, please post an issue if you get stuck or find a bug.

Thank you! I will try it when I have some time.