Carga crítica en una columna (G 7A)

clear all
L=9;
R0=1;
RL=5;
h=0.01;
N=(RL-R0)/h;
R=R0:h:RL;
P=zeros(1,N+1);
for k=1:(N+1)
    P(1,k)=(pi^3*1.2*R(1,k).^4)/2*L;
end
plot(R,P)
xlabel('Radio (m)');
ylabel('Carga crítica (KN)');

% Calculamos numéricamente los valores de P_cr. El primer autovalor es el que se corresponde con dicha carga (esto está almacenado en la primera columna de la matriz A)
clear all clf
a=0;
L=9;
ya=0;
yL=0;
h=0.1;
N=(L-a)/h;
x=a:h:L;
K=(1/h^(2))*(2*diag(ones(N-1,1))-diag(ones(N-2,1),1)-diag(ones(N-2,1),-1));
R=linspace(1,5,N+1);
for i=1:length(R)
    A(i,:)=eig(0.25*pi*1.2*R(i)^4.*K);
end
hold on
Pcr=A(:,1);
L2=9; p2=1.2;
R2=1:0.1:5;
P2=(pi^3*p2*R2.^4)/(4*L2.^2);
plot(R2,P2,'g')
plot(R,Pcr,'-')
xlabel('Radio (m)');
ylabel('Carga critica (kN)')
hold off