A fourth-order time-splitting Laguerre-Hermite pseudospectral method is introduced for Bose-Einstein condensates (BEC) in 3-D with cylindrical symmetry. The method is explicit, time reversible and time transverse invariant. It conserves the position density, and is spectral accurate in space and fourth-order accurate in time. Moreover, the new method has two other important advantages: (i) it reduces a 3-D problem with cylindrical symmetry to an effective 2-D problem; (ii) it solves the problem in the whole space instead of in a truncated artificial computational domain. The method is applied to vector Gross-Pitaevskii equations (VGPEs) for multi-component BECs. Extensive numerical tests are presented for 1-D GPE, 2-D GPE with radial symmetry, 3-D GPE with cylindrical symmetry as well as 3-D VGPEs for two-component BECs to show the efficiency and accuracy of the new numerical method.