The multitooth meshing state of harmonic drive (HD) is an important basic characteristic of its high transformation precision and high bearing capacity. Meshing force distribution affects the load sharing of the tooth during meshing, and theoretical research remains insufficient at present. To calculate the spatial distributed meshing forces and loading backlashes along the axial direction, an iterative algorithm and finite element model (FEM) is proposed to investigate the meshing state under varied transmission loading. The displacement formulae of meshing point under tangential force are derived according to the torsion of the flexspline cylinder and the bending of the tooth. Based on the relationship of meshing forces and circumferential displacements, meshing forces and loading backlashes in three cross-sections are calculated with the algorithm under gradually increased rotation angles of circular spline, and the results are compared with FEM. Owing to the taper deformation of the cup-shaped flexspline, the smallest initial backlash and the earliest meshing point appear in the front cross-section far from the cup bottom, and then the teeth in the middle cross-section of the tooth rim enter the meshing and carry most of the loading. Theoretical and numerical research show that the flexibility is quite different for varied meshing points and tangential force amplitude because of the change of contact status between the flexspline and the wave generator. The meshing forces and torsional stiffness of the HD are nonlinear with the torsional angle.