kill(all);

define unknown reaction forces at constraints, 
based on equilibrium

eq_tz  : WA+WC+WD-F  =0;
eq_rDx : (WA+WC)*b    =0;
eq_rDy : (-WC+F)*2*a=0;

solve the equilibrium equations

linsolve( [ eq_tz,eq_rDx,eq_rDy ] , [WD,WC,WA] ), globalsolve=true;

report the bending and the torsional moments for each beam segment;
see the attached figure for nomenclature and orientations

Mf1_AB : Mf1_CB : Mf1_DE : Mf1_FE : F*s$ Mf1_BE:0    $
Mf2_AB : Mf2_CB : Mf2_DE : Mf2_FE : 0  $ Mf2_BE:0    $
Mt__AB : Mt__CB : Mt__DE : Mt__FE : 0  $ Mt__BE:-2*a*F$

U:
 integrate( Mf1_AB^2/2/E/J + Mf2_AB^2/2/E/J +Mt__AB^2/2/G/Kt , s , 0 , a)
+integrate( Mf1_CB^2/2/E/J + Mf2_CB^2/2/E/J +Mt__CB^2/2/G/Kt , s , 0 , a)
+integrate( Mf1_DE^2/2/E/J + Mf2_DE^2/2/E/J +Mt__DE^2/2/G/Kt , s , 0 , a)
+integrate( Mf1_FE^2/2/E/J + Mf2_FE^2/2/E/J +Mt__FE^2/2/G/Kt , s , 0 , a)
+integrate( Mf1_BE^2/2/E/J + Mf2_BE^2/2/E/J +Mt__BE^2/2/G/Kt , s , 0 , b);

evaluate according to the second Castigliano theorem the deflection 
at the F force application point

d_F : diff(U,F);

evaluate the deflection according to the exercise numerical data

dim:    [a=800,b=1000,c=600,F=1000,C=500000,
	 E=210000,G=210000/2/(1+3/10),
	 J=(40^4-36^4)*%pi/64,Kt=(40^4-36^4)*%pi/32];

deflection, in mm

ev(d_F,dim,numer);