x=[1 2 3 4 5];
y=[0.1267e5 0.2219e5 0.2771e5 5.6971e5 2.2572e5];
plot(x,y,'-*')
xlabel('mode number'),ylabel('\omega^2_p')
title('P=\rho+3log\rho')