Numerical schemes for solving 3-D paraxial equations are constructed using splitting techniques. The solution can be reduced to a series of 2-D paraxial equations in each direction of splitting. The discretization along the depth is based on higher-order conservative schemes. The discretization along the transverse variables is based on higher-order finite difference variational schemes. Numerical experiments illustrate the advantages of higher-order schemes, which are much less dispersive, even for a small number of discretization points per wavelength.