Clear["Global`*"];

Beam bending using Ritz method
Consider a simply supported beam with a distributed load of the form: q(x) = q0 x1(L - x1)  and a concentrated load acting at x1 = (L/2). Compute an approximate solution. The properties of the beam are as follows:

1) Seek approximate solutions of the form:

Ck[n_] :=  Table[ck[i], {i, n}]
ϕk[n_]  :=  
      Table[Sin[(2 i  - 1) π (x1/L)], {i, 
          n}]
U3[n_]  :=  Ck[n] . ϕk[n]

U3[1]

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

d2U3[n_]  :=  
      D[U3[n], {x1, 2}]  // Simplify
U3qt[
        n_]  :=  U3[n]  /.  
        x1  ->  (L/4)
q  =  
        q0 x1 (L - x1);
Π[n_] :=(∫_0^L)( ((ym im)/2) (d2U3[n])^2 - 
                  U3[n] q) \[DifferentialD]x1 - P U3qt[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Π[1]

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

solution[2]

U3final[n_]  :=  U3[n]  /.  solution[n]

U3final[2]

Plot the solution

values  =  {ym  ->  70 (10^9), L ->  3, 
          im ->  ((0.05 (0.3^3))/12), q0 ->  -(10^5),  
          P ->   5 (10^5)};
Plot[Evaluate[U3final[3] /.  values], {x1, 0, L /.  values},  
        PlotRange -> All, 
        PlotStyle ->  {Thickness[0.01]}];

Beam deflection at x = (L/2)

values  =  {ym  ->  70 (10^9), L ->  3, 
          im ->  ((0.05 (0.3^3))/12), q0 ->  -(10^5),  
          P ->   4 (10^5)};
δ[
        n_]  :=  U3final[n]  /.  x1  ->  (L/2) /.  values // 
        N
convergence[
        n_]  :=  (δ[n + 1]  -  δ[n])/ δ[
          n + 1]
deflections   =  
      Table[δ[i],  {i, 
          4}]

Table[(((deflections[[n + 1]] - deflections[[n]]))/deflections)[[n]] , {n, 
        3}]