In this work we define, analyze, and compare different numerical schemes that
can be used to study the ground state properties of Bose-Fermi systems, such as
mixtures of different atomic species under external forces or self-bound
quantum droplets. The bosonic atoms are assumed to be condensed and are
described by the generalized Gross-Pitaevskii equation. The fermionic atoms, on
the other hand, are treated individually, and each atom is associated with a
wave function whose evolution follows the Hartree-Fock equation. We solve such
a formulated set of equations using a variety of methods, including those based
on adiabatic switching of interactions and the imaginary time propagation
technique combined with the Gram-Schmidt orthonormalization or the
diagonalization of the Hamiltonian matrix. We show how different algorithms
compete at the numerical level by studying the mixture in the range of
parameters covering the formation of self-bound quantum Bose-Fermi droplets.