kill(   
    all   
);

l : a-h;

I declare two auxiliary functions for linear interpolation along a segment

note: the evalution of the RHS of := is suspended (not performed) at define/assigment time; the suspended evaluation is performed in addition to the usual ones at call time

N1(xi) := 1-xi;

the define function, instead, evaluates as usual its second argument, then it assigns it to the first (which is NEVER evaluated, since object of further assigment, this is an exception of the general rule)

define(
    N2(goofy),
    goofy
);

define(
    mylinearinterpolationonBC(s),
    factor( vBp*N1(s/l)+vCm*N2(s/l) )
);

I initialize the strain energy as a null term

U:0$

Components of internal action along the AB segment, in-plane components alone

M1Ap:0*P+0*YA+0*XA-1*PsiA+0*F $

M1Bm:0*P+0*YA+h*XA-1*PsiA+0*F $

define(
    M1AB(s),
    expand(M1Ap*N1(s/h)+M1Bm*N2(s/h))
);

U : U + integrate( M1AB(s)^2/2/E/J , s , 0 , h );

BC segment

M1Bp:M1Bm $

M1Cm:0*P+0*YA+a*XA-1*PsiA-l*F $

define(
    M1BC(s),
    expand(M1Bp*N1(s/l)+M1Cm*N2(s/l))
);

U : U + integrate( M1BC(s)^2/2/E/J , s , 0 , l );

CD segment

M1Cp:M1Cm $

M1Dm:0*P-b*YA+a*XA-1*PsiA-l*F $

define(
    M1CD(s),
    expand(M1Cp*N1(s/b)+M1Dm*N2(s/b))
);

U : U + integrate( M1CD(s)^2/2/E/J , s , 0 , b );

DE segment

M1Dp:M1Dm $

M1Em:+a*P-b*YA+0*XA-1*PsiA+h*F $

define(
    M1DE(s),
    expand(M1Dp*N1(s/a)+M1Em*N2(s/a))
);

U : U + integrate( M1DE(s)^2/2/E/J , s , 0 , a );

U$
ratsimp(%)$
factor(%);

I proceed evaluating the partial derivatives of the strain energy with respect to the applied gen. loads in order to obtain the associated gen. displacements

diff(U,F)$
ratsimp(%)$
utB:%$

diff(U,P)$
ratsimp(%)$
utC:%$

diff(U,XA)$
ratsimp(%)$
uA:%$

diff(U,YA)$
ratsimp(%)$
vA:%$

diff(U,PsiA)$
ratsimp(%)$
psiA:%$

After its usage within the Castigliano theorem application, F may be finally declared as "fictitious"

F:0;

I now enforce the constraints at A by imposing null deflections there.

eqns : [ uA=0 , vA=0, psiA=0 ]$ 

unks : [ XA, YA, PsiA ];

implicit_ev( linsolve ( eqns , unks ) , globalsolve=true )

linsolve ( eqns , unks ) , globalsolve=true$

variant without "globalsolve"
sols:linsolve ( eqns , unks );
[XA,YA,PsiA]:ev([XA,YA,PsiA],sols);

it's time to substitute the obtained values for the constraint reactions at A within the sampled displacements at B and C.

ev([utB,utC])$
[utB,utC]:factor(%);

check that utC is not unduely dependent on h

diff(utC,h);
factor(%);

check that utB -> utC when h -> a

ratsimp( 
    limit( utB/utC , h , a )
);

I try to express the deflections in terms of a reference length (e.g. a) and some dimensional ratios ( r=b/a , t=h/a ) 

[utB,utC]$
ratsubst( t , h/a , % )$
ratsubst( r , b/a , % )$
factor(%)$
[utB,utC] : %;

utB;

Evaluation of a reference deflection obtained for a cantilever beam structure equivalent to a single upright frame member. 

bending moment for the ref. structure

Mfref:P*s$

strain energy for the ref. structure

Uref:integrate(Mfref^2/2/E/J,s,0,a);

I evaluate the deflection by resorting to the II Castigliano theorem

dref : diff(Uref,P);

I compare the frame stiffness (k=P/utC) with the stiffness of the one of its upright members alone, thus obtaining a "stiffness ratio"

ipkratio : (P/utC) / (P/dref);

I plot this normalized stiffness ratio

wxplot2d( 
    ipkratio , 
    [ r , 0.01 , 100 ] , 
    [logx] , 
    [legend, "ipkratio" ],
    [xlabel, "r=b/a ratio"],
    [ylabel, "ipkratio"]
);

ipkratio, r=1,numer;