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[(2 i  - 1) π (x/L)], {i, 1, 
          n}]
(Subscript[W, 0])[n_]  :=  Ck[n] . ϕk[n]

(Subscript[W, 0])[1]

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
(Subscript[W, tip])[n_]  :=  
      (Subscript[W, 0])[n]  /.  x  ->  L
Π[n_] :=(∫_0^L)( (1/2) ym im (d2Wdx[n])^2  +  
                  (NN/2) (((dWdx[n]))^2)) \[DifferentialD]x -  
          P (Subscript[W, tip])[n] // 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 length

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

matrix[2]

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

solution[n_]  :=  
      Solve[gradΠ[n]  ==  0 ,  Ck[n]]

solution[1]

(Subscript[W, approx])[n_]  :=  (Subscript[W, 0])[n]  /.  solution[n]

(Subscript[W, approx])[2]

Plot the solution

values  =  {ym  ->  70 (10^9), L ->  3, 
          im ->  ((0.05 (0.3^3))/12),  
          P ->   10000};
(Subscript[w, tip])[
        NNN_]  :=  (Subscript[W, approx])[2]   /.  x -> L  /.  
          values  /.  NN  ->  NNN
Plot[Evaluate[(Subscript[w, tip])[NNN]], {NNN, 0, -1950000},  
        PlotStyle ->  {Thickness[0.01]}];

(Subscript[W, approx])[3]   /.  
        x -> L /. values