Summary: The analog equation method, a BEM-based method, is employed for nonlinear analysis of Bernoulli-Euler beam with variable stiffness undergoing large deflections, under general boundary conditions which may be nonlinear. As the cross-sectional properties of the beam vary along its axis, the coefficients of the differential equations governing the equilibrium of the beam are variable. The formulation is in terms of displacements. The governing equations are derived in both deformed and undeformed configurations, and the deviations of the two approaches are studied. Using the concept of analog equation, the two coupled nonlinear differential equations with variable coefficients are replaced by two uncoupled linear ones pertaining to the axial and transverse deformation of a substitute beam with unit axial and bending stiffness, respectively, under fictitious load distributions. Besides the effectiveness and accuracy of the developed method, a significant advantage is that the displacements as well as the stress resultants are computed at any cross-section of the beam using the respective integral representations as mathematical formulae. Several beams are analyzed under various boundary conditions and loadings to illustrate the merits of the method as well as its applicability, efficiency and accuracy.