del[f_] := {D[f,bx],D[f,by],D[f,bz]} equilibrium = { bx -> 0 , by -> 0, bz -> 0} Vext = a(E^(2 bx) + E^(2 by) + E^(2 bz) ) Tad = b E^(-2 (bx+by+bz)/3 ) U = (xi -1 ) Tad ForceA = D[ Tad + U +Vext, bx] /. equilibrium bvalue = Simplify[ Solve[ForceA == 0 , b]][[1]] Ead = Simplify[ Tad +U +Vext /. bvalue ] KA = del[del[ Ead ]] /. equilibrium ihalf = { {1,0,0},{0,1,0},{0,0,lam}} mat = ihalf.KA.ihalf / ( a 2) omega2A = Eigenvalues[mat] Tdia = b (E^(-2 bx) + E^(-2 by) +E^( -2 bz) ) /3 Edia = Tdia + U + Vext /. bvalue Kdia = del[ del [ Edia ]] /. equilibrium mat = Simplify[ ihalf.Kdia.ihalf / ( 2 a) ] omega2dia = Eigenvalues[mat]