Fix boundary conditions

[timcallow]

The matrix diagonalization routine which is used to compute the KS orbitals and their eigenvalues does not give exactly the boundary conditions desired.

• For the dirichlet condition, it fixes the value of the eigenfunctions to zero at the ngrid +1 grid point; we want them equal to zero at the ngrid grid point
• For the neumann condition, it correctly sets the derivative to zero at the ngrid point, but is unable to resolve the eigenfunction at that point and returns a junk value

Although these only produced very small numerical errors, it is better to fix them. Their behaviour might also be more critical when we consider higher-order properties and more advanced features.

This PR introduces several fixes to ensure the correct behaviour of both bcs. Notably, the dirichlet grid is modified to first ensure the zero occurs at the ngrid point. For both bcs, the final point is determined via manual integration of the numerov formula.