A backward computation method has been developed to accelerate modelling of the critical state magnetization current in a staggered-array bulk high-temperature superconducting (HTS) undulator. The key concept is as follows: i) a large magnetization current is first generated on the surface of the HTS bulks after rapid field-cooling (FC) magnetization; ii) the magnetization current then relaxes inwards step-by-step obeying the critical state model; iii) after tens of backward iterations the magnetization current reaches a steady state. The simulation results show excellent agreement with the H-formulation method for both the electromagnetic and electromagnetic-mechanical coupled analyses, but with significantly faster computation speed. The simulation results using the backward computation method are further validated by the recent experimental results of a five period Gd-Ba-Cu-O (GdBCO) bulk undulator. Solving the finite element analysis (FEA) model with 1.8 million degrees of freedom (DOFs), the backward computation method takes less than 1.4 hours, an order of magnitude or higher faster than other state-of-the-art numerical methods. Finally, the models are used to investigate the influence of the mechanical stress on the distribution of the critical state magnetization current and the undulator field along the central axis.