A spectral collocation approach based on integrated polynomials is presented to investigate the statics and free vibrations of Euler–Bernoulli beams with axially variable cross section, modulus of elasticity, and mass density. The basic concept of the approach is the expansion of the highest derivatives appearing in the governing equations instead of the solution function itself by the truncated basis function. Then lower order derivatives and the function itself are obtained by integration. The constants appearing from the integrating process are determined by given classical or elastic restrained boundary conditions. Also, by incorporating the decomposition technique into the present approach, higher order vibration modes can be achieved even for stepped beams. Numerical examples including the statics and free vibrations of the beams with variance in geometry or material have been successfully solved, and the results are compared with those analytical or numerical solutions in the existing literature. The convergence and comparison studies show that convergent speed is rather rapid and the present approach can yield high accurate results with low computational efforts. Furthermore, the accuracy is not particularly affected by the adopted polynomials.