A search is performed for neutral non-standard-model Higgs bosons decaying to two muons in the context of the minimal supersymmetric standard model (MSSM). Proton-proton collision data recorded by the CMS experiment at the CERN Large Hadron Collider at a center-of-mass energy of 13 TeV were used, corresponding to an integrated luminosity of 35.9 fb$^{-1}$. The search is sensitive to neutral Higgs bosons produced via the gluon fusion process or in association with a $\mathrm{b\overline{b}}$ quark pair. No significant deviations from the standard model expectation are observed. Upper limits at 95% confidence level are set in the context of the $m_\mathrm{h}^{\text{mod+}}$ and phenomenological MSSM scenarios on the parameter $\tan\beta$ as a function of the mass of the pseudoscalar A boson, in the range from 130 to 600 GeV. The results are also used to set a model-independent limit on the product of the branching fraction for the decay into a muon pair and the cross section for the production of a scalar neutral boson, either via gluon fusion, or in association with b quarks, in the mass range from 130 to 1000 GeV.