You don't need to do a full FEM implementation for this problem. There are two simpler numerical approximations you can use. They are somewhat long and involved to write here, but I'll give you the references and a bird's eye view.
Use the beam-column frame theory. This is described in section 2.3 of Timoshenko and Gere's Theory of Elastic Instability. It requires you to get the beam properties and integrate the deflections and internal stresses. It's cumbersome, but not impossible.
Use the energy method. This is described in section 2.8 (and 2.9) of Timoshenko and Gere's Theory of Elastic Instability. In this method you basically "assume" a deformed shape (with some considerations the book treats) and then calculate the bending energy and the axial deformation energy. Buckling occurs when these energies are equal.
The simplified energy method is very simple, gives surprisingly accurate results, but requires heuristics to choose the approximated deformed shape.
I have personally implemented a partial version of 2.9 using Hermite polynomials of degree 3 (basically a simplified FEM formulation). I say partially because I only treated the simply supported case.
The advantage of this method is that, once implemented it will treat all frame cases and can be made as accurate as desired, by simply increasing the number of divisions of the element. The disadvantage is that it requires the calculation of the first eigenvalue of a matrix, which can be tricky if you don't have access to a linear algebra package.
In any case, read the 3 sections I mentioned (it's about 10-15 pages in total IIRC) and you can choose whatever method fits your needs best.
If you choose to go with the Hermite polynomial route let me know and I'll find some reference to that.