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

The finite element method V
Formulation of finite element matrices for an Euler-Bernoulli beam

Hermite Cubic Polynomials (used as basis functions)

(Subscript[ϕ, 1])  =  
        1  -  3 ((((ξ/h)))^2) +  
          2 ((((ξ/h)))^3);
(Subscript[ϕ, 2]) = ξ (((1  -  (ξ/h)))^2);
(Subscript[ϕ, 3]) =  
        3 ((((ξ/h)))^2) -  
          2 ((((ξ/h)))^3);
(Subscript[ϕ, 4])  =  ξ (((((ξ/h)))^2) -  (ξ/h));
ϕ  =  Table[(Subscript[ϕ, k]), {k, 4}];

ϕ  /.  ξ -> 0
(∂_ξ) ϕ  /.  ξ -> 
        0
ϕ  /.  ξ -> h
(∂_ξ) ϕ  /.  ξ -> h

Plot[Evaluate[ϕ /. h -> 1], {ξ, 0, 1},PlotStyle -> {{RGBColor[1, 0, 0], 
                  Thickness[0.01]}, {RGBColor[0, 1, 0], 
                  Thickness[0.01]}, {RGBColor[0, 0, 1], 
                  Thickness[0.01]}, {RGBColor[0, 1, 1], Thickness[0.01]}}];

Element stiffness matrix (Subscript[K, e])

Table[
      d2ϕ(Subscript[dx2, i]) = (∂_(ξ, ξ))(Subscript[ϕ, i]), {i, 
        4}]

(Subscript[K, e]) = 
      Table[(∫_0^h)( 
                EE II d2ϕ(Subscript[dx2, i]) 
                  d2ϕ(Subscript[dx2, j]) ) \[DifferentialD]ξ , {i, 4}, {j, 
            4}] // Simplify

((Subscript[K, e])/(2 EE (II/(h^3)))) // MatrixForm

Stiffness matrix should be singular, independently of h, E, I

Det[(Subscript[K, e])]

Element force vector

PP =  Table[P[i], {i, 4}]; 
    (Subscript[R, e]) = -(∫_0^h) (Subscript[q, 0]) ϕ  \[DifferentialD]ξ   +  PP // Simplify // MatrixForm