% This program corresponds to Appendix 3.2.3 of the thesis.
% The lateral stiffness of a Euler-Bernoulli column subjected to a tensile
%  axial load, T, is solved
clc, clear, close all

% Solution to the governing differential equation
syms C1 C2 phi T D L Q tl Rl tu Ru x EI phi
y(x) = C1*exp(phi*x/L) + C2*exp(-phi*x/L) - tl*Rl/T + Q*x/T;
dy(x) = diff(y,x);

% Substitute boundary conditions and external moment
eqn1 = y(0) == 0;
eqn2 = y(L) == D;
eqn3 = dy(0) == tl;
eqn4 = dy(L) == tu;
eqn5 = tu*Ru + tl*Rl == Q*L - T*D;

% Solve system of equations for Q
[A,B] = equationsToMatrix([eqn1,eqn2,eqn3,eqn4,eqn5],[Q,tu,tl,C1,C2]);
X = linsolve(A,B)
Q_f = simplify(X(1))

% Simplify Q_f manually by hand and verify the solution
ru = 1/(1+3*EI/(Ru*L));
rl = 1/(1+3*EI/(Rl*L));
a1 = 3*(ru*(1-rl)+rl*(1-ru));
a2 = 9*ru*rl+(phi)^2*(1-ru)*(1-rl);
a3 = 18*rl*ru - a1*(phi)^2;
Q_sol = D*12*EI/L^3*phi^3/12*(a1*phi*(exp(2*phi)+1) + a2*(exp(2*phi)-1)) /...
        ( 36*ru*rl*exp(phi) - a3*(exp(2*phi)+1) + (a2-a1)*phi*(exp(2*phi)-1) );
zero = simplify(subs(Q_sol-Q_f,phi,sqrt(T/EI)*L)) % Should be zero
