how to calculate beam buckling in matlab

main.m
L = 2;              % length of beam in meters
I = 4e-4;           % area moment of inertia in m^4
E = 200e9;          % modulus of elasticity in Pa
nu = 0.3;           % Poisson's ratio
187 chars
5 lines

main.m
% boundary conditions
bc = 'ss';          % simply supported
% load
q = 100;            % uniformly distributed load in N/m
124 chars
5 lines

main.m
% stiffness matrix
K = [12*E*I/L^3      6*E*I/L^2      -12*E*I/L^3     6*E*I/L^2;
     6*E*I/L^2       4*E*I/L        -6*E*I/L^2      2*E*I/L;
     -12*E*I/L^3    -6*E*I/L^2      12*E*I/L^3     -6*E*I/L^2;
     6*E*I/L^2       2*E*I/L        -6*E*I/L^2      4*E*I/L];
% load vector
P = [-q*L/2; q*L^2/12; -q*L/2; -q*L^2/12];
325 chars
8 lines

main.m
% eigenvalue analysis
[V,D] = eig(K,P);               % eigenvalue analysis
[lambda,idx] = sort(diag(D));   % sort eigenvalues in ascending order
V = V(:,idx);                   % corresponding eigenvector matrix
Pcr = -lambda(1)                % critical buckling load
phi = V(:,1)                    % mode shape
315 chars
7 lines

main.m
% plot mode shape
x = linspace(0,L,100);
y = [0; phi(2); 0; phi(4)]*sin(pi/L*x)';
plot(x,y)
title('First mode shape')
xlabel('x')
ylabel('y')
142 chars
8 lines