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) $$
