kill(all);

dominant term of the

f:sqrt(1+a^2)-1;

form

taylor( f ,a,0,2);

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
        +
        1/2 * (E*J) * k(myu(a,b),myv(a,b))^2,
    x,0,l
    )
    +
    F*ev(myu(a,b),x=l)
);

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

preloaded condition

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

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

[
    diff(%,da),
    diff(%,db)
];

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]),eta=1;

determinant(%);

linsolve(%=0,F);

ev(eqnl,%),eta=1;

linsolve(%,[da,db]);