restart; assume(a>0);
psi:=N*r*exp(-a*r);
                            psi := N r exp(-a r)
IN:=int(int(int(psi^2*r^2*sin(t),r=0..infinity),t=0..Pi),f=0..2*Pi);
                                         2
                                      3 N  Pi
                                IN := -------
                                         5
                                        a
N1:=solve(IN=1,N);
                   (1/2)       (1/2)  2     (1/2)       (1/2)  2
                  3      (Pi a)      a     3      (Pi a)      a
            N1 := ---------------------, - ---------------------
                          3 Pi                     3 Pi
N:=N1[1];
                               (1/2)       (1/2)  2
                              3      (Pi a)      a
                         N := ---------------------
                                      3 Pi
evalf(int(int(int(psi^2*r^2*sin(t),r=1/a..4/a),t=0..Pi),f=0..2*Pi));
                                0.8477145819
int(int(int(psi^2*r^1*sin(t),r=0..infinity),t=0..Pi),f=0..2*Pi); # 1/r
                                      a
                                      -
                                      2
int(int(int(psi^2*r^4*sin(t),r=0..infinity),t=0..Pi),f=0..2*Pi); # r^2
                                     15
                                    ----
                                       2
                                    2 a
int(int(int(psi*r^2*sin(t)*r^2*(sin(t))^2*psi,r=0..infinity),t=0..Pi),f=0..2*Pi);# x^2+y^2
                                     5
                                     --
                                      2
                                     a
int(int(int(psi*r^2*sin(t)*r*sin(t)*psi,r=0..infinity),t=0..Pi),f=0..2*Pi);# sqrt(x^2+y^2)
                                    5 Pi
                                    ----
                                    8 a