The first observation of the $D^0 to pi^+pi^-mu^+mu^-$ and $D^0 to K^+K^-mu^+mu^-$ decays is reported using a sample of proton-proton collisions collected by LHCb at a center-of-mass energy of 8$,$TeV, and corresponding to 2$,$fb$^{-1}$ of integrated luminosity. The corresponding branching fractions are measured using as normalization the decay $D^0 to K^- pi^+[mu^+mu^-]_{rho^0/omega}$, where the two muons are consistent with coming from the decay of a $rho^0$ or $omega$ meson. The results are $mathcal{B}(D^0 to pi^+pi^-mu^+mu^-)=(9.64pm0.48pm0.51pm0.97)times10^{-7}$ and $mathcal{B}(D^0 to K^+K^-mu^+mu^-)=( 1.54pm0.27pm0.09pm0.16)times10^{-7}$, where the uncertainties are statistical, systematic, and due to the limited knowledge of the normalization branching fraction. The dependence of the branching fraction on the dimuon mass is also investigated.