The drive-level dependence (DLD) of resonant frequency and the associated path-dependent, multi-valued amplitude-frequency response are critical nonlinear phenomena limiting the power handling and stability of high-frequency piezoelectric resonators. While Tiersten's seminal analytical model provides physical insight, it relies on an a priori introduction of an empirical quality factor (Q) to predict finite amplitudes, limiting its predictive power for design [1]. Subsequent finite-element analyses, though based on rigorous 3D nonlinear theory, face challenges in efficiently and robustly tracing the complete solution path, including unstable branches, across hysteretic regions [2]. This work presents a novel computational framework that bridges this gap. We develop a simplified, yet physics-based, 1D nonlinear boundary-value problem for the thickness-shear vibrations of monoclinic crystal plates (e.g., AT-cut quartz), incorporating viscous damping directly through material constants rather than an empirical Q [3]. Crucially, we implement a weighted pseudo-arclength continuation algorithm to solve this system. This method enables the automatic and continuous tracking of the full, path-dependent frequency-response curve under varying drive levels, seamlessly capturing jump phenomena and multi-valued regions without a priori assumptions. Our results demonstrate excellent agreement with established models for low drive levels while successfully predicting the complete hysteretic response at high drives. The framework operates solely on intrinsic material parameters (high-order elastic/piezoelectric coefficients, viscosity), requiring no device-specific fitting factors. It thus provides a more general, predictive, and computationally efficient tool for analyzing nonlinear resonance and stability in piezoelectric devices, surpassing the limitations of both prior analytical and numerical approaches.