  We develop a high accuracy power series method for solving partial
differential equations with emphasis on the nonlinear Schr\"odinger equations.
The accuracy and computing speed can be systematically and arbitrarily
increased to orders of magnitude larger than those of other methods. Machine
precision accuracy can be easily reached and sustained for long evolution times
within rather short computing time. In-depth analysis and characterisation for
all sources of error are performed by comparing the numerical solutions with
the exact analytical ones. Exact and approximate boundary conditions are
considered and shown to minimise errors for solutions with finite background.
The method is extended to cases with external potentials and coupled nonlinear
Schr\"odinger equations.
