The steady-state heat equation example from the MethodOfLines.jl documentation seems to break for smaller step sizes. The example I’m referring to is here: Steady State Heat Equation - No Time Dependence - NonlinearProblem · MethodOfLines.jl

Output for `dx == dy == 0.1`

(as in the original tutorial):

Output for `dx == dy == 0.09`

:

Output for `dx == dy == 0.05`

(all `NaN`

s on the interior):

The return code in every `PDESolution`

is `MaxIters`

(including the `0.1`

case), so I also tried significantly increasing `maxiters`

from the default. This doesn’t seem to do much:

Is there some reason why this is expected behavior?

Minimal example (essentially identical to the example in the docs, except that I’m using Makie for plotting and varying `dx`

/`dy`

and `maxiters`

):

```
using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve
using CairoMakie
@parameters x y
@variables u(..)
Dxx = Differential(x)^2
Dyy = Differential(y)^2
eq = Dxx(u(x, y)) + Dyy(u(x, y)) ~ 0
bcs = [u(0, y) ~ x * y,
u(1, y) ~ x * y,
u(x, 0) ~ x * y,
u(x, 1) ~ x * y]
# Space and time domains
domains = [x ∈ Interval(0.0, 1.0),
y ∈ Interval(0.0, 1.0)]
@named pdesys = PDESystem([eq], bcs, domains, [x, y], [u(x, y)])
dx, dy = 0.1, 0.1 # works great
#dx, dy = 0.09, 0.09 # instability near (1, 1)
#dx, dy = 0.08, 0.08 # unstable region grows
#dx, dy = 0.05, 0.05 # all NaN on the interior of the domain
iters_power = 5
maxiters = Int(10^iters_power)
# Note that we pass in `nothing` for the time variable `t` here since we
# are creating a stationary problem without a dependence on time, only space.
discretization = MOLFiniteDifference([x => dx, y => dy], nothing, approx_order=2)
prob = discretize(pdesys, discretization)
sol = NonlinearSolve.solve(prob, NewtonRaphson(), maxiters=maxiters)
begin
fig = Figure(resolution=(1000, 800))
ax = Axis(fig[1, 1], title="dx = $dx, dy = $dy, maxiters = 1e$iters_power")
h = heatmap!(ax, sol[x], sol[y], sol[u(x, y)], colorrange=(0, 1))
cb = Colorbar(fig[1, 2], h, label="u(x, y)")
fig
end
save("heat_ss_$(dx)_1e$(iters_power).png", fig)
```

Thanks in advance for taking a look!

Edit: Updated with more examples.