%
  
  function psi = energy_fun(P,q);
  
  alpha_1 = q(1,1);
  alpha_11 = q(2,1);
  alpha_111 = q(3,1);

  psi = alpha_1*P.^2 + alpha_11*P.^4 + alpha_111*P.^6;