  We present in details a numerical approach for solving supersymmetric quantum
mechanical systems with a gauge symmetry valid in all fermionic sectors. The
method uses a recursive algorithm to calculate matrix elements of any gauge
invariant operator in the Fock basis, in particular of the Hamiltonian
operator, and can be used for any gauge group. We describe its application to a
supersymmetric anharmonic oscillator model with discrete spectrum.
