Implementing such a method is not a trivial task. Just FYI, there are multiple programs purpose-built to perform Hartree-Fock calculations.
In terms of a good introduction to the theory of Hartree-Fock calculations, I found this pdf extremely helpful.
First of all, you have given the expression for the psuedo-exact (psuedo, because we have made the Born-Oppenheimer approximation and assumed a Slater determinant form to get here) 1-electron molecular orbitals $|\chi_a\rangle$, which are not known.
Formally, that's not a problem. We can just use some complete set of states $\{|\phi_\mu\rangle\}_{\nu=1...\infty}$ as a basis, and expand the exact (unknown) states by $|\chi_a\rangle = \sum_{\mu=1}^\infty C_{\mu a}|\phi_k\rangle$ in unknown coefficients $C_{\mu a}$.
Naturally, the requirement of infinitely many basis orbitals is not possible on a computer, so we must instead use a sufficiently large basis that the truncation doesn't matter to within our accuracy guidelines.
Following the analysis in the pdf, there are ultimately 4 matrices that we need.
$F_{\mu \nu}(C_{\alpha k}) = \langle \phi_\mu | \hat{f} | \phi_\nu \rangle$, the Fock matrix (The brackets are included here to represent dependence on $C$!)
$S_{\mu \nu}= \langle \phi_\mu | \phi_\nu \rangle $, the overlap integral of all of your basis functions (This does not change throughout the iteration.)
$C_{\alpha k}$, the coefficients that specify the expansion of the $k$th eigenfunction.
$\mathbf{\varepsilon}$, a diagonal matrix with entries corresponding to the 1-electron eigen-energies of $|\chi_a\rangle$ under Fock operator $\hat{f}$.
The choice of a basis "reduces" the earlier equation to
$$ \mathbf{F}(\mathbf{C}) \mathbf{C} = \mathbf{S} \mathbf{C} \mathbf{\varepsilon} $$
This type of equation is known as a generalised eigenvalue equation, and can be solved by most numerical linear algebra packages (e.g. Eigen, SciPy). Or rather, it would be if $F$ did not depend on $C$.
There are many, many algorithms for dealing with this, but a common approach is the "iteration" many authors refer to. In this setup, we start off with a guess for the form of $C$, get a better value of C and re-substitute and solve again.
$$ \mathbf{F}(\mathbf{C}^{(n)}) \mathbf{C}^{(n+1)} = \mathbf{S} \mathbf{C}^{(n+1)} \mathbf{\varepsilon} $$
where the superscripts on the $C$ matrices refer to the iteration number.
Why would the method necessarily converge
In general, it doesn't. We essentially hope that our starting guess was sufficiently close to the true solution that the algorithm can find the minimum. This is why it is of critical importance that the basis set is carefully chosen to be 'fairly close' to the true wavefunction from the get-go - A natural first choice is atomic hydrogenic wavefunctions, but these are computationally expensive to integrate and better results can be achieved with a larger collection of simpler functions. Such basis sets are tabulated for computational chemists at sites such as basis set exchange.
What sort of numerical method would I use to solve the Hartree-Fock equation?
Ultimately, it comes down to following this procedure:
- Choose a set of basis orbitals.
- Compute $\mathbf{S}$.
- Establish a guess for $C^{(0)}$.
- Solve $ \mathbf{F}(\mathbf{C}^{(n)}) \mathbf{C}^{(n+1)} = \mathbf{S} \mathbf{C}^{(n+1)} \mathbf{\varepsilon} $
- Repeat 4. until convergence.
A note on that convergence - generally, one tests to see if the new energies differ from the preceding energies by less than some cutoff, and so establish that the algorithm isn't "moving" very much in state space.
As for why it works at all - Conceivably, for sufficiently restricted conditions the step 4) above may define a contraction mapping $\mathcal{F} : M_{N\times M}(\mathbb{C}) \to M_{N\times M}(\mathbb{C}), \mathbf{C} \in M_{N\times M}(\mathbb{C})$, and so have convergence from the Banach fixed point theorem. For a general basis though, these requirements are not satisfied.
