# Boundary Condition Error in MethodOfLines

**URL:** <https://discourse.julialang.org/t/boundary-condition-error-in-methodoflines/137819>\
**Category:** Modelling & Simulations\
**Tags:** question, pde, methodoflines\
**Created:** [June 26, 2026, 8:14pm UTC](https://discourse.julialang.org/t/boundary-condition-error-in-methodoflines/137819 "2026-06-26T20:14:08Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![ducksoverip](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/ducksoverip/32/31967_2.png) [@ducksoverip](https://discourse.julialang.org/u/ducksoverip)\
**Post date:** [June 26, 2026, 8:14pm UTC](https://discourse.julialang.org/t/boundary-condition-error-in-methodoflines/137819/1 "2026-06-26T20:14:09Z")

</div>

I’m trying to use MethodOfLines to solve the time-independent Schrodinger equation for hydrogen, as a prelude to more complicated systems. The equation is as follows:

\nabla^2 \psi(\vec{r}) - \frac{1}{r}\psi(\vec{r}) = E\psi(\vec{r})

This problem has analytically known solutions and eigenvalues. Following the tutorials for the [steady-state heat equation](https://docs.sciml.ai/MethodOfLines/stable/tutorials/heatss/) and [1D Schrodinger equation](https://docs.sciml.ai/MethodOfLines/stable/tutorials/schroedinger/) in MethodOfLines.jl, I’ve written the following:

```julia-auto
using ModelingToolkit, MethodOfLines, OrdinaryDiffEq, DomainSets, NonlinearSolve

# Define parameters (indep. vars)
@parameters x, y, z, E 

# Define variables (known and unknown functions of params)
@variables psi(..)

# Define operators
Dxx = Differential(x)^2
Dyy = Differential(y)^2
Dzz = Differential(z)^2

V(x, y, z) = -1/sqrt(x^2 + y^2 + z^2)

# Define PDE
eq = [E*psi(x, y, z) ~ Dxx(psi(x, y, z)) + Dyy(psi(x, y, z)) + Dzz(psi(x, y, z)) + V(x, y, z)*psi(x, y, z)]

# set domain ranges and construct domain
xmin, ymin, zmin = 0, 0, 0
xmax, ymax, zmax = 10, 10, 10 

domains = [x in Interval(xmin, xmax), y in Interval(ymin, ymax), z in Interval(zmin, zmax)]

# set boundary conditions
bcs = [psi(0, 0, 0) ~ 0,
       psi(xmax, ymax, zmax) ~ 0]

# construct system, discretize
@named sys = PDESystem(eq, bcs, domains, [x, y, z], [psi(x,y,z)])
disc = MOLFiniteDifference([x => 100, y => 100, z => 100])
prob = discretize(sys, disc)

# solve
sol = NonlinearSolve.solve(prob, NewtonRaphson())

println(sol[E])

```

(Note: the operators and variables are specified in Cartesian coordinates, since the [documentation](https://docs.sciml.ai/MethodOfLines/stable/) says its restricted to Cartesian grids. The same page says it can handle spherical Laplacians but then says nothing else and provides no examples, so go figure.) As written, this currently errors on the boundary condition `psi(0, 0, 0) ~ 0`, with the following message:

```julia-auto
ERROR: AssertionError: Boundary condition psi(0, 0, 0) ~ 0 is not on a boundary of the domain, or is not a valid boundary condition
Stacktrace:
 [1] parse_bcs(bcs::Vector{…}, v::PDEBase.VariableMap, orders::Dict{…})
   @ PDEBase C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\parse_boundaries.jl:370
 [2] symbolic_discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; checks::Bool)
   @ PDEBase C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\symbolic_discretize.jl:27
 [3] symbolic_discretize
   @ C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\symbolic_discretize.jl:9 [inlined]
 [4] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…}; analytic::Nothing, checks::Bool, kwargs::@Kwargs{})
   @ PDEBase C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\discretization_state.jl:74
 [5] discretize(pdesys::PDESystem, discretization::MOLFiniteDifference{…})
   @ PDEBase C:\Users\njsjr\.julia\packages\PDEBase\mILTI\src\discretization_state.jl:69
 [6] top-level scope
   @ c:\Users\njsjr\Documents\attosecond-research\helium-tdse\_research\pde_test.jl:36
Some type information was truncated. Use `show(err)` to see complete types.

```

I’m very confused by this, because I set the boundaries of my domain to be 0 and 10, so I don’t see why there should be an issue here.

---

<div class="post-metadata">

**Author:** ![technocrat](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/technocrat/32/220947_2.png) [@technocrat](https://discourse.julialang.org/u/technocrat)\
**Post date:** [June 26, 2026, 9:23pm UTC](https://discourse.julialang.org/t/boundary-condition-error-in-methodoflines/137819/2 "2026-06-26T21:23:10Z")

</div>

Bear with me. I hope to be useful even if it turns out that I’m misguided.

The immediate problem is that the `MethodOfLines` requires boundary conditions on faces, not vertices. `psi(0,0,0) ~ 0` defines a point in 3D, not a surface.

But, once past that **E** is an unknown Eigenvector and as used you’ll just get `ψ = 0` which is not very informative. What I _think_ is needed is to take an initial whack at a trial value for **E** or add a normalization constraint and treat **E** as a free variable but be careful to avoid singularity at the origin by a slight shift. in one dimension, say `xmin=0.1.` Take care not to start with a large grid size or you may find yourself wondering if it will ever finish. Better yet explolit spherical symmetry and reduce to one dimension in r,

```julia
-u''(r) + V(r)·u(r) = E·u(r)

```

which avoids the 3D grid.

But please don’t follow this rabbit too far into Wonderland, because I haven’t tested it. Offered just as a potentially useful insight.

---

<div class="post-metadata">

**Author:** ![tobydriscoll](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tobydriscoll/32/1843_2.png) [@tobydriscoll](https://discourse.julialang.org/u/tobydriscoll)\
**Post date:** [June 29, 2026, 12:53pm UTC](https://discourse.julialang.org/t/boundary-condition-error-in-methodoflines/137819/3 "2026-06-29T12:53:59Z")

</div>

What you have is an eigenvalue problem, which requires a different numerical approach than the two examples you followed. I’m not aware that MOL can handle eigenvalue problems.

For 1-D problems (including spherical symmetry, which is _not_ the case for a box), you could look at [ApproxFun](https://juliaapproximation.github.io/ApproxFun.jl/stable/generated/Eigenvalue/).