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. The search is sensitive to neutral Higgs bosons produced via the gluon fusion process or in association with a 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 and phenomenological MSSM scenarios on the parameter 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.