We investigate the prospects for the discovery of neutral Higgs bosons with muons by direct searches at the CERN Large Hadron Collider (LHC) as well as by indirect searches in the rare decay Bs -> mu+ mu- at the Fermilab Tevatron and the LHC. Promising results have been found for the minimal supersymmetric standard model, the minimal supergravity (mSUGRA) model, and supergravity models with non-universal Higgs masses (NUHM SUGRA). For tanb simeq 50, we find that (i) the contours for a branching fraction of B(Bs -> mu+ mu-) = 1x10^{-8} in the parameter space are very close to the 5sigma contours for pp -> b phi^0 -> b mu+ mu- + X, phi^0 = h^0, H^0, A^0 at the LHC with an integrated luminosity (L) of 30 fb^{-1},(ii) the regions covered by B(Bs -> mu+ mu-) ge 5x10^{-9} and the discovery region for bphi^0 -> b mu+ mu- with 300 fb^{-1} are complementary in the mSUGRA parameter space,(iii) in NUHM SUGRA models, a discovery of B(Bs -> mu+ mu-) simeq 5x10^{-9} at the LHC will cover regions of the parameter space beyond the direct search for pp -> bphi^0 -> b mu+ mu- with L = 300 fb^{-1}.