Abstract: The existing algorithms for Prolate Spheroidal Wave Functions (PSWFs) have poor efficiency, high complexity in hardware implementation and especially uncontrollable precision. To overcome the above weaknesses, a new algorithm based on state transition matrix approximation of differential equation is proposed combining the theory of linear time-varying system. In the new algorithm, the state transition matrix on the whole interval is approximated by the ones on very small intervals. After that, the movement track of the system on discrete time spot is attained, and that is the numerical solution of PSWFs. The expression of error to the precise value is deduced theoretically and then the algorithm is improved to get the briefer error expression. The new algorithm is compared and analyzed with the one proposed by Parr and the approximation algorithm of Legendre polynomials on the calculation precision and complexity. The simulation results show that the proposed algorithm has high and controllable precision, low complexity and are easy for hardware implementation.