eigenval.m

function [f1,f2,f3] = eigenval(N,P)
E = 70*10^9; %modulous of elasticity
d = 2700; %density
b = 0.1; %base
h = 0.01; %height
L = 1; %length
I = (b*(h^3))/12;
Area = b*h; 
segl = L/N;
Diag = ((P*2)/((N^2)*E*I)) + 6; %Diag values
OFFDiag = (-P/((N^2)*E*I)) - 4; 
A = zeros(N-1,N-1); % creates A matrix
for row = 1:(N-1)
    for column = 1:(N-1) %diag values
        if row == column
            A(row,column) = Diag;
        end
        if column == row - 1 %off diag values
            A(row,column) = OFFDiag;
        end
        if column == row + 1
            A(row,column) = OFFDiag;
        end
        if column == row - 2   %"off-off" diag values
            A(row,column) = 1;
        end
        if column == row + 2
            A(row,column) = 1;
        end
    end
end
A(1,1) = ( (2*P) / ((N^2)*(E*I)) ) + 5;%corner values set one less than inner values
A(N-1,N-1) = ( (2*P) / ((N^2)*(E*I)) ) + 5;
eigenval = eig(A);
dx = L/(N+1);
x = linspace(dx,L-dx,N+1);%requires natural freqs.
f1 = sqrt((eigenval(1)*E*I)/d/Area/(segl^4));
f2 = sqrt((eigenval(2)*E*I)/d/Area/(segl^4));
f3 = sqrt((eigenval(3)*E*I)/d/Area/(segl^4));
setdefaults
figure(1)
plot(x,sin(pi*x),x,sin(2*pi*x),x,sin(3*pi*x))
xlabel('x distance, m')
ylabel('Deflection')
end