Clear["Global`*"];
Off[General::spell, General::spell1];

Beam-column/buckling using Ritz method

1) Seek approximate solutions of the form:

Ck[n_] :=  Table[ck[i], {i, n}]
ϕk[n_]  :=  
      Table[Sin[i π (x/L)], {i, 1, n}]
(Subscript[W, 0])[n_]  :=  Ck[n] . ϕk[n]

(Subscript[W, 0])[3]

Plot[Evaluate[ϕk[4]  /.  L -> 1],  {x, 0, 1.2}]

2) Define the approximate potential Π(Ck) for this problem

dWdx[n_]  :=  
      D[(Subscript[W, 0])[n], x]  // Simplify
d2Wdx[
        n_]  :=  D[(Subscript[W, 0])[n], {x, 2}]  // Simplify
Π[n_] :=(∫_0^L)( (1/2) EI (d2Wdx[n])^2  +  
                (NN/2) (((dWdx[n]))^2)) \[DifferentialD]x // 
        Expand

Π[1]

3) Minimize the potential to find the value of the coefficients that best approximate this problem

gradΠ[n_]  :=  
      Table[D[Π[n], ck[i]], {i, n}]

gradΠ[2]

Buckling load

matrix[n_]  := 
      Table[Coefficient[gradΠ[n][[i]], Ck[n]], {i, 
          n}]

matrix[2]

Solve[Det[matrix[2]] == 0,  NN]

The second result is the critical load corresponding to the first buckling mode.