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.4 $\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.