A search for neutral non-standard-model Higgs bosons decaying to two muons is presented. The search is performed in the context of the minimal supersymmetric standard model, using proton-proton collisions data recorded by CMS at the CERN Large Hadron Collider at a center-of-mass energy of 13~\mathrm{TeV}. The integrated luminosity is 35.9~\mathrm{fb}^{-1}. The search is sensitive to neutral Higgs bosons produced via gluon fusion process or in association with a quark pair. No significant deviation from the standard model expectation is observed. A 95 confidence level upper limit is set in the context of the and hMSSM scenarios on the parameter as a function of the pseudoscalar A boson mass, in the range from to 600~\mathrm{GeV}. The larger collected luminosity and the higher center-of-mass energy exclude a larger - region, compared to what was obtained at and 8~\mathrm{TeV} by a similar analysis. 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 to 1000~\mathrm{GeV}.