General Solution to Linear Schrodinger equation I am trying to find a solution to
$$\displaystyle \left[-\frac{1}{2}\nabla^2 - \frac{2}{r} + C(r)\right]\phi(r) = E\phi(r)$$
where $C(r)$ is a known function of r. I am just looking for some help on the strategy to solve explicitly for $\phi(r)$. I have read through the solution for hydrogen using separation of variables but
1) I don't really understand it, they make a lot of handwavy arguments about terms disappearing  :)
2) How does one handle the extra $C(r)$?
FYI: this is for the Hartree Fock method on Helium so the actual equation is the one posted in the answer here: Hartree Fock equations
P.S. I posted this in math.stackechange first should I close that one? 
 A: I feel like this question is a little too open ended and I'm sorry for that...but I did find a WONDERFUL HF walk through http://www.phys.sinica.edu.tw/TIGP-NANO/Course/2011_Spring/classnotes/CMS_20110511.pdf that clears everything up.
A: To do Hartree-Fock to Helium atom, you just need to calculate one orbital, which for helium is spherically symmetric. The Hartree-Fock 'integro-differential' equation for spherically symmetric atom with one eigenstate can be written as
$$\left(-\frac{1}{2}\nabla^2 - \frac{2}{r} + V_{Hx}(r) \right) u(r) = \epsilon u(r),$$
where $u(r) = \psi(r) r$ and 
$$V_{Hx}(r) = \frac{1}{2} \int d{\bf r}' \frac{n({\bf r'})}{|{\bf r-r}'|}$$
The factor $\frac{1}{2}$ comes from exchange term cancelling half of the Hartree-potential.
I would recommend shooting method. (Alternatives: gaussian basis set, finite element, finite difference methods). 
Shooting method is simple. Guess an eigenvalue and integrate the solution of radial Schrödinger equation from r=0 to r=r_max. Require that boundary condition is fulfilled at r_max (zero). If not, change the eigenvalue accordingly.
https://en.wikipedia.org/wiki/Shooting_method
Solving the radial Hartree-Fock potential equation is also super simple. 
$$ \int dr' \frac{n(r')}{|r-r'|} = \frac{1}{r} \int_0^r dr' 4\pi n(r) r^2 + \int_r^\infty 4\pi n(r) r $$
(Eq. 7.6 here) https://wiki.fysik.dtu.dk/gpaw/_static/rostgaard_master.pdf
The latter part of my answer here explains how this equation works.
Dipole in a spherical cavity in an infinite dielectric
Beware, that the Fock-operator is only multiplicative for the Helium 1s orbitals. Rest of the eigenvalue spectrum of the given equation does not match to Fock-operator spectrum.
edit: I just read a comment of yours, which stated that you are using gaussian orbitals. You can solve it easily. You just need to calculate some matrix elements. Given a set of gaussian orbitals $\phi_i$, you calculate two matrices
$$H_{ij} = < \phi_i | H | \phi_j>$$ 
and 
$$S_{ij} = < \phi_i | \phi_j >$$.
Then you solve this generalized eigenvalue equation
$$Hc=\epsilon Sc$$,
and you get the new coefficients for the wave function as the eigenvector. Your wave function is now
$$ \psi(r) = \sum_n c_n \phi_n(r) $$
The method is called subspace diagonalization. 
To solve the generalized eigenvalue equation with octave or Matlab,
octave:1> H = [ -2 0.1 ; 0.1 -3 ];
octave:2> S = [ 1 0.1 ; 0.1 1 ];
octave:5> [psi, lambda] = eig(H,S)
psi =

  -0.35088  -0.94180
   0.97217  -0.25494

lambda =

Diagonal Matrix

  -3.1498        0
        0  -1.9209

Here you obtain the eigenvalues (-3.14 is the lowest), and the corresponding eigenvector [-0.35, 0.97]. Now, all that there is left, is to fill matrices H and S with the real matrix elements.
To be more spesific, following integrals must be evaluated:
$$H_{ij} = \int dr \phi_i(r) (-1/2 \nabla^2 + V(r) ) \phi_j(r) \\
S_{ij} =  \int dr \phi_i(r) \phi_j(r) $$
