I am currently developing a numerical solver, which solves the standard, one-dimensional, time-independent Schrödinger equation
$$\frac{-\hbar^2}{2m}\frac{d^2\psi}{dx^2} + V(x)\psi = E\psi \tag{1}$$ I want to make it such that my solver works for non-uniform grids. More precisely, I would like to be able to have one grid spacing $a_{1}$ for say $x<0$ and another $a_{2}$ for $x\geq 0$. To do this I employ the following finite difference discretization of (1)
$$\frac{-\hbar^2}{2m} \Bigl( \frac{\psi_{i+1}-\psi_{i}}{x_{i+1}-x_{i}} - \frac{\psi_{i}-\psi_{i-1}}{x_{i}-x_{i-1}} \Bigr)\cdot \frac{2}{x_{i+1}-x_{i-1}} + V_{i}\psi_{i} = E\psi_{i}\tag{2}$$
This can then be solved as an eigenvalue equation $\hat{H}\boldsymbol{\psi} = E\boldsymbol{\psi}$, and it turns out that this delivers consistent results, also for non-uniform grid spacings.
There is however one thing that makes me uncertain about the validity of (2). Namely that when applying it to a non-uniform grid the resulting Hamiltonian matrix $\hat{H}$ is not hermitian, which is evident if one considers a point lying on the boundary between grid points with spacing $a_{1}$ and grid points with spacing $a_{2}$. In this case $\hat{H}$ will have a entries, which do not match their corresponding hermitian conjugate entries.
Is this a problem? And if so, what would be the correct way to discretize (1) and preserve hermiticity?
