Simple Matlab Program to calculate Euler Buckling Load using Finite Element Method

This program determine the Euler buckling loads for simply supported column

EI d^2 y/ dx^2 +Py =0, y(0)=0, y(L)=0.
A weak form using weighing function w(x) such that w(0)=0, w(L)=0. w(x) is piece wise linear interpolation function.
Writing the weak form, we get
-int 0 to L dw/dx dy/dx dx + lamda int 0 to L w(x) y(x) dx =0
Using the linear shape functions as shown below:
w(x) = N1(x) w1+ N2(x) w2,
y(x)= N1(x) y1_+ N2(x) y2.
we get
kdet= k0*P/(E*I)-k1;
We find the values of P such that the det(kdet) goes to zero.