x=[1 2 3 4 5];
y=[5.4513e5 2.6741e5 -0.9480e5 -0.6233e5 -0.6233e5];

plot(x,y,'-*')
xlabel('mode number'),ylabel('\omega^2_p')
title('P=\rho+8log\rho')