Boundary conditions in axes?

EDIT: This was my original motivation, see next post for a more generalized take on the problem.

The Schrödinger equation in spherical coordinates

is singular at the origin. The solution has to obey the following boundary conditions:

In finite-differences, it is common to adjust the upper-left corner of the derivative matrices for \ell = 0 to implement this and preserve the convergence order (see refs below). I wonder how I could accomplish this with ContinuumArrays.jl, namely I wish to dispatch an expression like

R = ...
D = Derivative(axes(R,1))
Δ = R'D'D*R

correctly such that I get the correct boundary conditions? One idea could be to introduce

abstract type AbstractDerivative{T} <: LazyQuasiMatrix{T} end

struct Derivative{T,D} <: AbstractDerivative{T}
    axis::Inclusion{T,D}
end
...

# Then I could have my own type
struct SingularDerivative{T} <: AbstractDerivative{T}
    axis::Inclusion{T,D}
end

• If R is aware of SingularDerivative
  • If \ell == 0, apply correction
  • Else, return uncorrected matrix
• If R is not aware of SingularDerivative, proceed as normal

The problem with this approach is that suddenly all basis libraries will need to dispatch on AbstractDerivative.

A better alternative would of course be if I could somehow attach boundary conditions to the derivative operator:

D = Derivative(axes(R,1), :singular, :regular)

Are there better approaches?

References

• Schafer, K. J., Gaarde, M. B., Kulander, K. C., Sheehy, B., &
  DiMauro, L. F. (2000). Calculations of strong field multiphoton
  processes in alkali metal atoms. AIP Conference Proceedings, 525(1),
  45–58. http://dx.doi.org/10.1063/1.1291925

• Schafer, K. J. (2009). Numerical methods in strong field physics. In
  T. Brabec (Eds.), (pp. 111–145). : Springer.

• Muller, H. G. (1999). An Efficient Propagation Scheme for the
  Time-Dependent Schrödinger equation in the Velocity Gauge. Laser
  Physics, 9(1), 138–148.

• Patchkovskii, S., & Muller, H. (2016). Simple, accurate, and
  efficient implementation of 1-electron atomic time-dependent
  Schrödinger equation in spherical coordinates. Computer Physics
  Communications, 199(nil),
  153–169. http://dx.doi.org/10.1016/j.cpc.2015.10.014

Start by reading the definitions of Derivative, LazyQuasiMatrix, Inclusion, and the axes(R, 1) usage described in the issue, then inspect how basis libraries currently dispatch on derivative operators. The issue does not name files or tests; done would require an agreed design and documented behavior for singular boundary conditions without forcing every basis library to add bespoke dispatch.