egv.m

function [f1,f2,f3] = egv(N,P)
%Function solves for the natural frequency of a the beam in problem 2 and 3
%given N segments and transverse load P

%Initialize Constants
E = 70*10^9; %Pa
dens = 2700; %kg/m^3
b = 0.1; %m
h = 0.01; %m
l = 1; %m

%Derived variables
I = (b*(h^3))/12; % m^4
area = b*h; %m^2

%Segments
segl = l/N; %m - length of segment

%Diagonal values
diagonal = ((P*2)/((N^2)*E*I)) + 6; %main diagonal
offDiag = (-P/((N^2)*E*I)) - 4; %off-diagonal

%Loop through all values of A to create values
A = zeros(N-1,N-1);
for row = 1:(N-1)
    for column = 1:(N-1)
        %diagonal values
        if row == column
            A(row,column) = diagonal;
        end
        %off diagonal values
        if column == row - 1
            A(row,column) = offDiag;
        end
        if column == row + 1
            A(row,column) = offDiag;
        end
        %"off-off" diagonal values
        if column == row - 2
            A(row,column) = 1;
        end
        if column == row + 2
            A(row,column) = 1;
        end
    end
end

%Sets corner values to one less than main diagonal values
A(1,1) = ( (2*P) / ((N^2)*(E*I)) ) + 5;
A(N-1,N-1) = ( (2*P) / ((N^2)*(E*I)) ) + 5;

%Built in eigenvalue MATLAB function finds eigenvalues of the matrix
%generated for the beam
egv = eig(A);

%Plug into function for natural frequencies to solve (for first three -
%more can be generated with following values of egv vector)
f1 = sqrt((egv(1)*E*I)/dens/area/(segl^4));
f2 = sqrt((egv(2)*E*I)/dens/area/(segl^4));
f3 = sqrt((egv(3)*E*I)/dens/area/(segl^4));
end