© 2019 The Author(s). 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 13TeVwere used, corresponding to an integrated luminosity of 35.9fb-1. The search is sensitive to neutral Higgs bosons produced via the gluon fusion process or in association with a bbquark pair. No significant deviations from the standard model expectation are observed. Upper limits at 95% confidence level are set in the context of the mmod+hand phenomenological MSSM scenarios on the parameter tanβas a function of the mass of the pseudoscalar Aboson, in the range from 130 to 600GeV. 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 bquarks, in the mass range from 130 to 1000GeV.