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); /* -diff(v,x,2) should behave the same... */

e(u,v) and k(u,v) are the generalized strains obtained
from the displacements, i.e. they are the counterpart
of the 
[e k]^T  = B(u,v)*[u v]^T (???)

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, i.e. the sum of
- the elastic strain energy
- the unexerted work of external forces; 
  if assumed =0 in undeformed configuration, it is defined as
  the "Le" performed work of the external forces in reaching the 
  configuration under scrutiny, changed in sign (i.e. =-Le).

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)

[ maybe, lambda*F0 could be put in place of F, for added consistency
with LPBA ]

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.
in the geometrically nonlinear case

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

I simplify the equations by considering small perturbations

taylor(%,[da,db],[0,0],1);
ratdisrep(%); /* the "..." placeholder for the further terms is removed*/
eqnl:%$

the equations are homogeneus, as it is expected since
the unperturbed state (da=0,db=0) is an equilibrium condition

digression on [aug]coefmatrix maxima functions:

A*x=b
A*x-b=0

augcoefmatrix(A*x-b,[x])
[A -b]

   coefmatrix(%,[x])
[A]

In order to have nontrivial solution, 
the coefficient matrix of the eqnl equations
must be singular

coefmatrix(eqnl,[da,db]);

determinant(%);

I search for the "F" load which zeroes the determinant

defFcrit : linsolve(%=0,F);

and I search the eqnl solutions for that specific "F" value

%pi^2,numer;

eqnlcrit : ev(eqnl,defFcrit);

linsolve(eqnlcrit,[da,db]);

thus obtaining an arbitrary growth 
of the (supposedly small..) "db" perturbation
in the absence of applied load variation.