This looks promising, there is the time dependent Schroedinger equation (2) and the first order differential Eqs. (7) to (9) which can be solved simultaneously by Maxima. The numerical results can be checked with the analytical results. The radial wavefunctions psi(r) are the modified Laguerre polynomials and the angular psi (theta, phi) are the spherical harmonics. This looks like a new ab initio method which can be straightforwadly extended to helium without using perturbation theory. In the relativistic case it seems that the equivalents of Eqs. (7) to (9) are also soluble with Maxima for hydrogen. Then spin can be introduced with the Pauli matrices, and at a final stage the Dirac type wavefunctions introduced. The great advantage is the use of the Euler Lagrange equations to give r dot, theta dot and phi dot.