The Born-Oppenheimer (BO) approximation is less accurate in the presence of a strong magnetic field than in the absence of a field. This is due to the complicated and unpredictable response of electronic structure to the field, especially in the mixed regime $B \approx B_0 = 2.35 \times 10^5\,$T. Therefore, it is desirable to explore non-BO methods in magnetic fields. In this work, the nuclear-electronic orbital (NEO) method is employed to study protonic vibrational excitation energies in the presence of a strong magnetic field. NEO Generalized Hartree-Fock theory and time-dependent Hartree-Fock theory are derived and implemented, accounting for all terms that result as a consequence of the nonperturbative treatment of molecular systems in a magnetic field. The NEO results for HCN and FHF$^-$ with clamped heavy nuclei are compared against the quadratic eigenvalue problem (QEP). Each molecule has three semi-classical modes owing to the hydrogen - two precession modes that are degenerate in the absence of a field and one stretching mode. The NEO-TDHF model is found to perform well - in particular it automatically captures the screening effects of the electrons on the nuclei, which are quantified through the difference in energy of the precession modes.