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

% Solve the governing differential equation
syms EI Rl tl Q P kAG
syms y(x) psi(x)
ode1 = EI*diff(psi) == Rl*tl - Q*x + P*y
ode2 = Q - P*diff(y) == kAG*(-psi + diff(y)) % replace diff(y) with psi as necessary
odes = [ode1; ode2];
S = dsolve(odes);
y = simplify(S.y)
psi = simplify(S.psi)

% Verify the solutions for y(x) and psi(x) 
syms D L tu Ru
syms C1 C2 phi w
y(x) = C1*exp(phi*w/L*x) + C2*exp(-phi*w/L*x) - Rl*tl/P + Q/P*x;
psi(x) = Q/P + 1/w*(C1*phi/L*exp(phi*w/L*x) - C2*phi/L*exp(-phi*w/L*x));

% Substitute boundary conditions
eqn1 = y(0) == 0;
eqn2 = y(L) == D;
eqn3 = psi(0) == tl;
eqn4 = psi(L) == tu;
eqn5 = tu*Ru + tl*Rl == -P*D + Q*L;

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

% Simplify Q 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+(w*phi)^2*(1-ru)*(1-rl);
a3 = 18*rl*ru - a1*(w*phi)^2;
wphi = phi*w;
beta = (w*phi^3)/12 * ( a1*wphi*(exp(2*wphi)+1) + a2*(exp(2*wphi)-1) ) / ( 36*rl*ru*exp(wphi) - a3*(exp(2*wphi)+1) + (a2-a1)*wphi*(exp(2*wphi)-1) );
zeta = (w^2-1)*(36*rl*ru*exp(wphi) - 18*rl*ru*(exp(2*wphi)+1) - a1*wphi*(exp(2*wphi)-1))/   ( 36*rl*ru*exp(wphi) - a3*(exp(2*wphi)+1) + (a2-a1)*wphi*(exp(2*wphi)-1) );
Q_result = 12*EI/L^3*beta/(1+zeta) * D;
Q = subs(Q,phi,sqrt(P/EI)*L);
Q_result = subs(Q_result,phi,sqrt(P/EI)*L);

% Compare results:
error = simplify(Q_result-Q) % should be zero