kill(all);

PRELIMINARY RESULT

assume(dx>0);

I consider two neighboring points P and Q along the beam axis
P: position x, disps (uP,vP)
Q: position x+dx disps (uQ,uQ) defined as

uQ:uP+dudx*dx$ vQ:vP+dvdx*dx$

their initial distance l0 equates

PQ0:dx$

whether their distance in the deformed configuration equates

PQ:sqrt((dx+uQ-uP)^2+(vQ-vP)^2);

the axial elongation is defined as

etmp : factor( PQ/PQ0 ) - 1;

whose taylor expansion at the SECOND order is

etmp : taylor(etmp,[dudx,dvdx],[0,0],2);

END OF THE PRELIMINARY RESULT

kill(all);

large rotation formulations for elongation and curvatures (eta=1); 
the customary small rotation is obtained when the second order terms are neglected
(eta=0).

e(u,v):= diff(u,x) + 1/2*diff(v,x)^2 *eta;

k(u,v):=diff(v,x,2);

2 dof representation of the structure deformation

myu(a,b):=-a*x/l+b*0;

myv(a,b):=a*0+b*sin(x/l*%pi);

alternatively, if the exact solution is not known
we may employ the polynomial form
myv(a,b):=4*b*(x/l)*(1-x/l);

assume(  l>0  );

overall potential energy

define(
    U(a,b),
    integrate(
        1/2 * (E*A) * e(myu(a,b),myv(a,b))^2 /* e=N/EA */
        +
        1/2 * (E*J) * k(myu(a,b),myv(a,b))^2,/* k=Mf/EJ */
    x,0,l
    )
    +
    F*ev(myu(a,b),x=l) /* = -F*a */
);

equilibrium configuration are searched for as stationary points of the potential
energy in (a,b)

eq:[
    diff(U(a,b),a),
    diff(U(a,b),b)
];

preloaded condition
(aka. solution of the problem under linearity assumptions)

ev(eq,eta=0);
linsolve(%,[a,b]);
[a0,b0]:ev([a,b],%)$

the small rotation solution is a solution also for the large rotation problem

[a0,b0];

eq, a=a0,b=b0;

I consider a perturbation in the neighborhood of such preloaded condition

U( a0 +da , b0 +db );

and I search for equilibrium conditions for this perturbed config.

[
    diff(%,da), 
    diff(%,db)
],eta=1;

I simplify the equations by considering small perturbations

taylor(%,[da,db],[0,0],1);
ratdisrep(%);
eqnl:%$

A*x=b
A*x-b=0
augcoefmatrix(%,[x])
[A -b]

coefmatrix(eqnl,[da,db]);

determinant(%);

defFcrit : linsolve(%=0,F);

eqnlcrit : ev(eqnl,defFcrit);

linsolve(eqnlcrit,[da,db]);