eigenval.asv

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; %Diagonal values
OFFDiag = (-P/((N^2)*E*I)) - 4; 
A = zeros(N-1,N-1);%Loop through values of A to generate values
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" diagonal 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;%Sets corner values to 1 less than main diag 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);%Plug in for natural frequencies to solve
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