You should execute this line and the next before you try using this routine. restart;assume(a>0); The following finds the radial wave functions for hydrogen. a is the Bohr radius, r is the radial coordinate, and n and l are the principal and total angular momenta quantum numbers. radial:=proc(n,l) add((-1)^(i-l)*(2*r/n/a)^i/(i-l)!/(i+l+1)!/(n-i-1)!,i=l..n-1)*exp(-r/a/n);simplify(%/sqrt(int(%^2*r^2,r=0..infinity))) end proc: For example, if you want the radial wave function for n = 3 and l = 1, you just plug in these two values radial(3,1); LCQqLEkickc2IiIiIiwmSSNhfGlyR0YlIiInRiQhIiJGJkYoIyEiKCIiIy1JJGV4cEc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkdGJTYjLCQqJkYkRiZGKEYqI0YqIiIkRiZGKSNGJkYtI0YtIiRWIw== Is it normalized? Let's check it. int(radial(3,1)^2*r^2,r=0..infinity); IiIi Now let's set up a routine for spherical harmonics. We'll let l be the total angular momentum quantum number and m be the magnetic quantum number, or Lz eigenvalue. spherharm:=proc(l,m) local i;sin(theta)^l:for i from -l to m-1 do simplify(diff(%,theta)-cot(theta)*i*%) end do:%*sqrt((2*l+1)*(l-m)!/4/Pi/(l+m)!)/2^l/l!*exp(I*m*phi) end proc: For example, to find the spherical harmonic for l = 7 and m=4, just plug these numbers in spherharm(7,4); LCQqLi1JJHNpbkc2JCUqcHJvdGVjdGVkR0koX3N5c2xpYkc2IjYjSSZ0aGV0YUdGKSIiJS1JJGNvc0dGJkYqIiIiLCYqJEYtIiIjIiM4ISIkRi9GLyIkcSgjRi9GMkkjUGlHRicjISIiRjItSSRleHBHRiY2IyomXiNGLEYvSSRwaGlHRilGL0YvIyIiJCIjaw== Yuck! Is this even normalized properly? assume(phi::real);2*Pi*integrate(abs(spherharm(7,4))^2*sin(theta),theta=0..Pi); IiIi The product of these pairs of functions is then the hydrogen wave functions. Maple sometimes has some odd ideas of how to express the answers, however radial(5,3)*spherharm(3,1); LCQqNkkickc2IiIiJCwmSSNhfGlyR0YlIiM/RiQhIiIiIiJGKCMhIzYiIiMtSSRleHBHNiQlKnByb3RlY3RlZEdJKF9zeXNsaWJHRiU2IywkKiZGJEYrRihGKiNGKiIiJkYrIiNxI0YrRi4sJiokLUkkY29zR0YxNiNJJnRoZXRhR0YlRi5GOEYqRitGKy1JJHNpbkdGMUY/RisiI0BGOkkjUGlHRjIjRipGLi1GMDYjKiZeI0YrRitJJXBoaXxpckdGJUYrRisjRioiKCtEYyc=