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 NaNs) 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.