  A new numerical treatment in the Crank-Nicholson method with the imaginary
time evolution operator is presented in order to solve the Schr\"{o}dinger
equation. The original time evolution technique is extended to a new operator
that provides a systematic way to calculate not only eigenvalues of ground
state but also of excited states. This new method systematically produces
eigenvalues with accuracies of eleven digits with the Cornell potential that
covers non-perturbative regime. An absolute error estimation technique based on
a power counting rule is implemented. This method is examined with exactly
solvable problems and produces the numerical accuracy down to 10$^{-11}$.
