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 bb‾ 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 mhmod+ and phenomenological MSSM scenarios on the parameter tan⁡β 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.