The magnetic order of plutonium hydrides (PuHx) has long been a subject of controversy, both experimentally and theoretically. The magnetic, structural, electronic, and thermodynamic properties of PuHx are investigated using Hubbard-corrected density functional theory (DFT + U), with U derived from linear response calculations. To address the issue of electronic metastable states within DFT + U, we employ an occupation matrix control method as well as allowing 5f orbitals to break the structural symmetry. This advanced computational approach establishes that the ground-state magnetic order for PuH2 is antiferromagnetic and that for PuH3 is ferromagnetic, consistent with Faraday and nuclear magnetic resonance experiments. Using the conventional hydrogen-vacancy model for PuHx, the hydrogen-content-induced magnetic transition and anomalous variation of the magnetic moment are successfully reproduced. In particular, the calculated transition point of x matches the experimental findings. Furthermore, the electronic configuration of the Pu atom in PuHx is determined to be 5f5, which is consistent with observations from X-ray photoemission spectroscopy. Our calculated values for enthalpy of formation, heat capacity, and entropy are in good agreement with experimental data, thereby validating the robustness of our theoretical framework.