```function [phi,wn] = eigenModes(geometry,BC,Nmodes)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                               GOAL

%  Compute eigen frequencies and mode shapes of a beam with different
%  configurations and geometries

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                               INPUT

% geometry: structure with fields defined in the Example.m file
% BC : 1,2,3 or 4 : Boundaries conditions (cf Example.m)
% Nmodes : Number of modes wished

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%                              OUPUT
%                         ----------------
% Phi : Matrix of mode shapes is a [Nmodes x Nz] vector.
% wn is a  vector of size [1 x Nmodes ] in rad/sec.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%  check input

if ~isnumeric(BC),
error('BC must be an integer between 1 and 4');
end

if max(ismember([1,2,3,4],BC))==0,
error('BC is unknown. Please, choose BC = 1,2,3 or 4');
end

if ~isnumeric(Nmodes),
error('Nmodes must be an integer');
end

name = [{'L'},{'E'},{'nu'},{'rho'},{'I'},{'y'},{'m'}];
for ii=1:numel(name),
if ~isfield(geometry,name{ii}),
error([' The field ',name{ii},' is missing in the structure "geometry".'])
end
end

%%
y= geometry.y;
Nz = numel(y); % number of discrete elements

% volume V and mass m of the beam
m = geometry.m;
L = geometry.L;

% get the non trivial solution of f
Ndummy=1:Nmodes^2; % the number is arbitrary fixed as the square of the number of Nmodes.
h =fsolve(@modeShape,Ndummy);
% small values of h come from numerical errors
h(h<0.1)=[];
% Analytically, many solutions are identical to each other, but numerically, it is not the
% case. Therefore I need to limit the prceision of the solutions to 1e-4.
h = round(h*1e4).*1e-4;
% the uniques solution are called beta:
beta = unique(h);
beta = beta(1:Nmodes);

%modes shapes calculations
if Nmodes>numel(beta),
error(' Mode shape required is too high \n\n')
end
wn = beta.^2.*sqrt(geometry.E*geometry.I./(m.*geometry.L^4)); % in rad
phi = zeros(Nmodes,Nz);

for ii=1:Nmodes,
switch BC,

case 1,% pinned-pinned
phi(ii,:)=sin(beta(ii)*y./L);

phi(ii,:)=phi(ii,:)./max(abs(phi(ii,:)));

case 2 % clamped-free

phi(ii,:)=(cosh(beta(ii)*y./L)-cos(beta(ii)*y./L))+...
((cosh(beta(ii))+cos(beta(ii))).*(sin(beta(ii)*y./L)-sinh(beta(ii)*y./L)))./...
(sin(beta(ii))+sinh(beta(ii)));

phi(ii,:)=phi(ii,:)./max(abs(phi(ii,:)));

case 3 % clamped-clamped

phi(ii,:)=sin(beta(ii)*y./L)-sinh(beta(ii)*y./L)-...
((sin(beta(ii))-sinh(beta(ii))).*(cos(beta(ii)*y./L)-cosh(beta(ii)*y./L)))./...
(cos(beta(ii))-cosh(beta(ii)));

phi(ii,:)=phi(ii,:)./max(abs(phi(ii,:)));

case 4 % clamped-pinned

phi(ii,:)=sin(beta(ii)*y./L)-sinh(beta(ii)*y./L)-...
((sin(beta(ii))-sinh(beta(ii))).*(cos(beta(ii)*y./L)-cosh(beta(ii)*y./L)))./...
(cos(beta(ii))-cosh(beta(ii)));

phi(ii,:)=phi(ii,:)./max(abs(phi(ii,:)));
end

end

% Function to solve is different for each BC:

function [f] = modeShape(y)
switch BC,
case 1,% pinned-pinned
f=sin(y).*sinh(y);
case 2 % clamped-free
f=cos(y).*cosh(y)+1;
case 3 % clamped-clamped
f=cos(y).*cosh(y)-1;
case 4 % clamped-pinned
f=cos(y).*sinh(y)-sin(y).*cosh(y);
end
end

end

```