From ab0922800ce2ec55c8bd2cc94bd3c8ad7aec7cf0 Mon Sep 17 00:00:00 2001
From: Shuhei Ohno <ohno.github@gmail.com>
Date: Fri, 14 Aug 2026 12:33:01 +0900
Subject: [PATCH] Apply parity boundary conditions

---
 docs/src/FDM.md | 10 ++++++++++
 src/FDM.jl      | 17 +++++++++++++++--
 test/FDM.jl     | 14 ++++++++++++++
 3 files changed, 39 insertions(+), 2 deletions(-)

diff --git a/docs/src/FDM.md b/docs/src/FDM.md
index 8f9f65f..7eff3bc 100644
--- a/docs/src/FDM.md
+++ b/docs/src/FDM.md
@@ -21,6 +21,16 @@ The eigenvalue ``E`` is an approximation of the exact energy and the eigenvector
 ```
 A uniform grid spacing is used, ``r_{i+1} = r_{i} + \Delta r``. See the API reference for the expression of the matrix ``\pmb{H}``.
 
+## Boundary conditions
+
+For central differences, the wavefunction is extended across the origin with
+
+```math
+\psi(-r)=(-1)^l\psi(r).
+```
+
+Thus even ``l`` uses ``\psi'(0)=0``, while odd ``l`` uses ``\psi(0)=0``. The origin value for even ``l`` is eliminated with a second-order one-sided derivative, so the radial grid can continue to exclude ``r=0``.
+
 ## Usage
 
 Run the following code before each use.
diff --git a/src/FDM.jl b/src/FDM.jl
index 6fe77c8..d522552 100644
--- a/src/FDM.jl
+++ b/src/FDM.jl
@@ -52,11 +52,24 @@ function matrix(o::RestEnergy, method::FiniteDifferenceMethod)
 end
 
 function matrix(o::Kinetic, method::FiniteDifferenceMethod)
-  D  = FiniteDifferenceMatrices.fdmatrix(Int64(length(method.R)), n=1, m=2, d=method.direction, h=method.Δr, t=typeof(method.Δr))
-  D² = FiniteDifferenceMatrices.fdmatrix(Int64(length(method.R)), n=2, m=2, d=method.direction, h=method.Δr, t=typeof(method.Δr))
+  D, D² = _derivative_matrices(method)
   return -o.hbar^2/2/o.m * (D² + SparseArrays.spdiagm(2 ./ method.R) * D - method.l*(method.l+1) * SparseArrays.spdiagm(1 ./ method.R .^ 2))
 end
 
+function _derivative_matrices(method::FiniteDifferenceMethod)
+  n = Int64(length(method.R))
+  D = FiniteDifferenceMatrices.fdmatrix(n, n=1, m=2, d=method.direction, h=method.Δr, t=typeof(method.Δr))
+  D² = FiniteDifferenceMatrices.fdmatrix(n, n=2, m=2, d=method.direction, h=method.Δr, t=typeof(method.Δr))
+
+  if method.direction == :c && iseven(method.l) && 1 < n
+    D[1,1] = -2 / (3method.Δr)
+    D[1,2] = 2 / (3method.Δr)
+    D²[1,1] = -2 / (3method.Δr^2)
+    D²[1,2] = 2 / (3method.Δr^2)
+  end
+  return D, D²
+end
+
 function matrix(o::PotentialTerm, method::FiniteDifferenceMethod)
   return SparseArrays.spdiagm([V(o, r) for r in method.R])
 end
diff --git a/test/FDM.jl b/test/FDM.jl
index 01b9834..2ec4f90 100644
--- a/test/FDM.jl
+++ b/test/FDM.jl
@@ -1,5 +1,19 @@
 @testset "FDM.jl" begin
 
+  @testset "boundary condition" begin
+    even = FiniteDifferenceMethod(Δr=0.1, rₘₐₓ=0.3, l=0)
+    odd = FiniteDifferenceMethod(Δr=0.1, rₘₐₓ=0.3, l=1)
+    D_even, D²_even = TwoBody._derivative_matrices(even)
+    D_odd, D²_odd = TwoBody._derivative_matrices(odd)
+
+    even_wavefunction = 1 .+ even.R .^ 2
+    odd_wavefunction = collect(odd.R)
+    @test (D_even * even_wavefunction)[1] ≈ 2even.R[1]
+    @test (D²_even * even_wavefunction)[1] ≈ 2
+    @test (D_odd * odd_wavefunction)[1] ≈ 1
+    @test (D²_odd * odd_wavefunction)[1] ≈ 0 atol=1e-12
+  end
+
   # Testing Results
 
   H = Hamiltonian(