# MethodOfLines - numerical instability for example in docs

**URL:** https://discourse.julialang.org/t/methodoflines-numerical-instability-for-example-in-docs/106359
**Category:** New to Julia
**Tags:** modelingtoolkit, nonlinearsolve, methodoflines
**Created:** [November 17, 2023, 2:58pm UTC](https://discourse.julialang.org/t/methodoflines-numerical-instability-for-example-in-docs/106359 "2023-11-17T14:58:20Z")
**Posts on this page:** 2
**Page:** 1

<div class="post-metadata">

### Author: ![Sushrut\_Deshpande](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sushrut_deshpande/32/52754_2.png) [@Sushrut\_Deshpande](https://discourse.julialang.org/u/Sushrut_Deshpande)
#### Post date: [November 17, 2023, 2:58pm UTC](https://discourse.julialang.org/t/methodoflines-numerical-instability-for-example-in-docs/106359/1 "2023-11-17T14:58:20Z")

</div>

I am using the steady state heat equation example in the MethodOfLines.jl package.

When I reduce the step side below 0.1 I have a numerical blow-up of the solution.

I do not know if this is only on my end or other users have also faced this error.  
To me I do not see a reason for numerical instability of this case. I have attached my code below:

```julia
using ModelingToolkit, MethodOfLines, DomainSets, NonlinearSolve

@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 = 0.05
dy = 0.05

# 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())

u_sol = sol[u(x, y)]

using Plots

heatmap(sol[x], sol[y], u_sol, xlabel="x values", ylabel="y values",
        title="Steady State Heat Equation")

```

This is the solution I get:  
 ![NumericalBlowup](https://global.discourse-cdn.com/julialang/original/3X/c/a/cabf3a310add5309c7e7c817235dd727ef7ce93a.png)

---

<div class="post-metadata">

### Author: ![xtalax](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/xtalax/32/35293_2.png) [@xtalax](https://discourse.julialang.org/u/xtalax)
#### Post date: [November 18, 2023, 5:57pm UTC](https://discourse.julialang.org/t/methodoflines-numerical-instability-for-example-in-docs/106359/2 "2023-11-18T17:57:57Z")

</div>

Yes, I have heard of this before. Difficult to say what might be causing it, there’s nothing special about 0.1 as a number