A search for decays of the Higgs and $Z$ boson to pairs of $J/\psi$ or $\Upsilon(nS) (n = 1, 2, 3)$ mesons, with their subsequent decay to $\mu^+$ $\mu^−$ pairs, is presented. The analysis uses data from proton-proton collisions at $\sqrt{s} = 13$ TeV, collected with the CMS detector at the LHC in 2017 and corresponding to an integrated luminosity of 37.5 fb$^{−1}$ . While an observation of such a decay with this sample would indicate the presence of physics beyond the standard model, no significant excess is observed. Upper limits at 95% confidence level are placed on the branching fractions of these decays. In the $J/\psi$ pair channel, the limits are $1.8 \times 10^{-3}$ and $2.2 \times 10^{-6}$ for the Higgs and Z boson, respectively, while in the combined $\Upsilon(nS)$ pair channel, the limits are $1.8 \times 10^{-3}$ and $1.5 \times 10^{-6}$ , respectively, when the mesons from the Higgs and Z boson decay are assumed to be unpolarized. When fully longitudinal and transverse polarizations are considered the limits reduce by about $22–29$% and increase by about $10–13$%, respectively.