Plotting basis functions

All
 I am trying to plot the 6-31G basis functions for carbon along one
 coordinate (eventually two for p functions).
 'Trying' is the right word...
 I am using the cartesian expression given in Levine, 5th, p 487, which
 appears to be the general form of those found in other texts. The basis
 set information comes from the EMSL library.
 About the only thing that works is that all the functions go to zero at
 large distance from the nucleus. Looking along one axis, the 1s (6G)
 function looks correct. However, the 2s (3G) function has a probability
 near zero at the nucleus while the 2p function does not have a zero
 probability at the nucleus!
 The 1s(6G), 2s(3G), and 2p(3G) _do not_ integrate to unity!?. They are
 close, between 0.86 and 1.13, but this is too large for errors in the
 integration routine (MathCad, from -5 a0 to 5 a0 along all coordinates).
 As a check, the 1G functions do integrate to unity.
 Any hints would be appreciated!
 I presume the 6-31G set was not optimized on hydrogenic orbitals, but
 shouldn't there be some resemblance?
 Roy Jensen
 -----BEGIN PGP PUBLIC KEY BLOCK-----
 Comment: available at www.consol.ca
 ------END PGP PUBLIC KEY BLOCK------