format long
% compute the five natural frequencies of a free-free beam
%
% ENDS YOUR INPUT DATA. THE REST DOES IT FOR YOU!
%
betaL_try = [0,4.6, 7.5,12, 14 ];
betaL = zeros(1,0);
for n=1: 5,
%  k = 10^n; % the support spring constants relative to the tension T, k_L = k_0 = k;
%
% compute the first natural frequency
%
x =betaL_try(n);  % the starting value from completely free-free case;
norm = 1;
while (norm>1.e-07)

A =  (1 - cos(x)*cosh(x));
norm = abs(A);

slope = sin(x)*cosh(x) - cos(x)*sinh(x);      

if (abs(slope)>1.e-08)
xnew = x - A/slope;
else
xnew = x;
end;
x= xnew;
end;
betaL=[betaL x];

% compute the first natural frequency in Hertz
%omega = (betaL^2/L^2)* sqrt(EI/m);
%freq = omega/(2*pi); %in Hertz
%
% compute the mode shape

% characteristic matrix, eq. 23 of lecture 12 note
% A  = [
%  0  -1  0  1;
% -1  0   1   0;
% -s  -c  sh  ch;
% -c   s  ch  sh];
%
s = sin(x); c=cos(x); sh=sinh(x); ch = cosh(x);
Aminor=[ 
  0   -1  0;
 -1   0   1;
 -s  -c  sh];

 b=[-1; 0; -ch];
coef = Aminor\b;


x_coord=zeros(1,0);
y_coord =zeros(1,0);

for i=0:100,
span = i/100;
arg =x*span;
sx = sin(arg); cx = cos(arg);
shx =sinh(arg); chx = cosh(arg);
W =coef(1)*sx + coef(2)*cx + coef(3)*shx + chx; 
x_coord = [x_coord span];
y_coord =[y_coord W];
end;

%normalize y_coord
y_coord = y_coord/(max(abs(y_coord)));
figure(n);
plot(x_coord, y_coord);
xlabel('Beam  span');
ylabel('Mode shape ');
legend([ num2str(n), '-th  mode ',  '  beta*L =  ', num2str(x)]);
title('Mode shape of free-free beam');

end;  %end the outer loop


disp('   ');
disp(['beta L value of first five modes  ', num2str(betaL)]);