A quick and convenient way to allow for the effect of deformations on the behaviour of beams and columns under axial load is to adjust the beam stiffness matrix, so that the calculated deflections are in equilibrium with the applied loads, including the eccentricity of the axial load due to the transverse deflections.  Further background to this approach is presented in the book linked in the previous post:  Matrix Structural Analysis.

The geometric stiffness matrix for a 1D beam (from the link above) is:

I have incorporated this approach into a spreadsheet that will calculate deflections, shear forces and bending moments for a beam subject to combined axial load and a uniform transverse load.  The spreadsheet, including full open source code, may be downloaded from:

Geom_K.zip

Note that the spreadsheet uses the Strand7 axes convention, which changes the sign of some of the matrix values:

Input and output for the Axbeam function is shown below.  The beam may be subdivided into any number of segments, each of which may have different materials properties. Deflections or shear forces and bending moments are output at the end of each segment.

Typical input and output is shown below for a beam with 16 segments, and 2 different material properties.  Results are shown for an axial load of 1000 kN, and also 1800 kN, which is close to the buckling load:

Output deflections are shown below for the two load cases, compared with results for the same beam in Strand7.  For the lower load the results are almost identical.  For the higher load the iterative procedure used in Strand7 has resulted in significantly smaller deflections:

A similar trend is seen for bending moments and shear forces:

The MaxDefAx function returns the maximum absolute deflection for a series of axial loads. In the example below the beam has the same EI value over the full length.  It can be seen that the predicted buckling load matches the Euler buckling value very closely:

The beam end conditions may be adjusted.  In the example below the moment fixity has been removed at both ends.  The spreadsheet buckling load again matched the Euler value almost exactly (the maximum value in the table has been reduced below the buckling load, to show the shape of the deflection curve more clearly).  When one or more end fixities are released, end loads may be applied to that freedom.