Problem 9.11 from textbook

The finite element model corresponds to two cubic (hermite) (C^1)interpolation beam elements.
       (Subscript[e, 1])          	(Subscript[e, 2])
o------------o------------o
\(TraditionalForm\`\(\(U\_1\) \(U\_2\)\(\ \ \ \ \ \ \)\(U\_3\) \(U\_\(\(4\)\(\ \ \ \ \ \ \ \)\)\)\(\ \ \ \ \ \ \ \)\(U\_5\)\(\ \)\(U\_6\)\(\ \ \)\)\)
where (Subscript[U, 2]),  (Subscript[U, 4]),  (Subscript[U, 6])are rotations. 

The  stiffness matrix of a generic element is:
(K^e)=((2  (E^e) (I^e))/(h^3))(({{6, -3 h, -6, -3 h}, {-3 h, 2 (h^2), 3 h, (h^2)}, {-6, 3 h, 6, 3 h}, {-3 h, (h^2), 3 h, 2 (h^2)}}))
The variationally-consistent load vector corresponding to a uniformly distributed load is:

R=(({{-((h (Subscript[q, 0]))/2)}, {-(((h^2) (Subscript[q, 0]))/12)}, {-((h (Subscript[q, 0]))/2)}, {(((h^2) (Subscript[q, 0]))/12)}}))

The global (assembled) stiffness matrix is obtained as follows:

h = 0.12; (Subscript[I, 1]) =  N[π ((0.03^4)/4), 30]; 
    (Subscript[E, s]) =  200 (10^9);
K1=N[((2  (Subscript[E, s]) (Subscript[I, 1]))/(h^3))(({{6, -3 h, -6, -3 h}, {-3 h, 2 (h^2), 3 h, (h^2)}, {-6, 3 h, 6, 3 h}, {-3 h, (h^2), 3 h, 2 (h^2)}})),30]
(Subscript[I, 2]) =  π ((0.02^4)/4);(Subscript[E, s]) =  200 (10^9);K2=((2  (Subscript[E, s]) (Subscript[I, 2]))/(h^3))(({{6, -3 h, -6, -3 h}, {-3 h, 2 (h^2), 3 h, (h^2)}, {-6, 3 h, 6, 3 h}, {-3 h, (h^2), 3 h, 2 (h^2)}}))

K = (({{K1[[1]][[1]], K1[[1]][[2]], K1[[1]][[3]], K1[[1]][[4]], 0, 0}, {K1[[2]][[1]], K1[[2]][[2]], K1[[2]][[3]], K1[[2]][[4]], 0, 0}, {K1[[3]][[1]], K1[[3]][[2]], K1[[3]][[3]] + K2[[1]][[1]], K1[[3]][[4]] + K2[[1]][[2]], K2[[1]][[3]], K2[[1]][[4]]}, {K1[[4]][[1]], K1[[4]][[2]], K1[[4]][[3]] + K2[[2]][[1]], K1[[4]][[4]] + K2[[2]][[2]], K2[[2]][[3]], K2[[2]][[4]]}, {0, 0, K2[[3]][[1]], K2[[3]][[2]], K2[[3]][[3]], K2[[3]][[4]]}, {0, 0, K2[[4]][[1]], K2[[4]][[2]], K2[[4]][[3]], K2[[4]][[4]]}}));
MatrixForm[
      K]

The force vector is:

h  =  0.12; (Subscript[q, 0]) =  200; (Subscript[P, 3]) =  1000;  
    (Subscript[M, 6]) =  5000;
r= (({{(Subscript[R, 1])}, {(Subscript[R, 2])}, {-h ((Subscript[q, 0])/2) - (Subscript[P, 3])}, {-(h^2) ((Subscript[q, 0])/12)}, {(Subscript[R, 5])}, {(Subscript[M, 6])}}))//Flatten;
r

The solution vector is defined as:

u = (({{0}, {0}, {(Subscript[U, 3])}, {(Subscript[U, 4])}, {0}, {(Subscript[U, 6])}}))//Flatten;

The solution is obtained doing:

Kreduced=(({{K[[3, 3]], K[[3, 4]], K[[3, 6]]}, {K[[4, 3]], K[[4, 4]], K[[4, 6]]}, {K[[6, 3]], K[[6, 4]], K[[6, 6]]}}))
Rreduced =  {r[[3]],  
        r[[4]],  r[[6]]}

Inverse[Kreduced] . Rreduced

Alternatively

solution  =  
      Solve[K . u == r,  {(Subscript[R, 1]), (Subscript[R, 2]),  (Subscript[U, 3]),  (Subscript[U, 4]),  (Subscript[R, 5]), (Subscript[U, 6])}] // 
        Flatten