Tesseral perturbation is a major part in analytical orbit propagation, which involves a significant amount of calculations of inclination functions and their partial derivatives. Recursion formula is widely applied for accuracy and efficiency purposes. Based on the modified Gooding's method proposed in the literature, an optimized program is proved to reduce 24% of calculation time for all inclination functions up to degree/order 50 ($L\le 50$). Furthermore, considering that the mean inclination varies slightly in orbit prediction for 1--3 days, the functions are expanded as Taylor series, which speeds up calculations significantly, and the time consumed during orbit prediction is reduced by 48% while the tesseral perturbation is calculated up to degree/order 50. Moreover, with expansion of inclination functions up to the second order, the predicted ephemerides deviate slightly from those calculated using full recursion formula. For a 500km low Earth satellite, propagating its orbits for 3 days using modified Gooding's method and second-order Taylor expansion respectively to calculate the inclination functions within degree/order 50, the ephemeris deviation RMS (Root Mean Square) is less than 1mm and decreases as the altitude increases.