kill(all);

I name "1" the supported node, and I replace the constraint with its reaction "R": subsequently I will impose the compatibility condition of null vertical (z) displacement at this node.
To define linear functions along segments, I use a normalized curvilinear abscissa xi, which
it is 0 at the first point of the segment, and 1 at the final point.

only a bending moment component is nonzero, being the structure planar, and loaded purely out-of-plane

Mf_12: -R*a*xi;
Mf_23: (F-R)*b*xi;

torsional moment

Mt_12: 0;
Mt_23: -R*a;

I integrate the elastic strain energy contributes along the various beam segments.

U: integrate( (Mf_12^2/2/E/J + Mt_12^2/2/G/Jp) *a, xi, 0, 1)
  +integrate( (Mf_23^2/2/E/J + Mt_23^2/2/G/Jp) *b, xi, 0, 1);

I employ the Castigliano theorem to evaluate the deflection at the application point of the "R" force

delta_R: diff(U,R);

I impose the compatibility condition described above, and I obtain the value of "R".

linsolve(delta_R=0,R),globalsolve=true;

I redefine U by substituting the expression of R obtained above, thus expliciting its dependence on F.

U: ev(U);

By skipping this apparently redundant passage, the final derivative  of U in F might erroneously neglect the contributions hidden within the unsubstituted R symbols.

I reduce to a common denominator form, to shorten the formula

U:fullratsimp(U);

I evaluate the deflection at the application point of F

delta_F: diff(U,F);

I define a specific set of dimensions, to derive a numerical evaluation of such a deflection

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,Jp=(40^4-36^4)*%pi/32];

ev(delta_F,dim)$
%,numer;

Such a result may be compared with the output of a consistent FE analysis.