J=1;
delta=1;
kappa=0.1*J;
eta=15*J; g0=10^-4*J; gamma=0.002*J; omega_m=0.1*J; K=0.05*J; h=6.626*10^-34;
hbar=h/(2*pi);
kB=1.38*10^-23;
t=0;
nbar=(exp(hbar*omega_m/(kB*t))-1)^-1;
syms bet alpha
format compact
eqns=[bet-g0*abs(alpha)^2/(omega_m-1i.*gamma-2*K)==0,alpha-1i*eta/(delta+1i*kappa+2*J+2*g0*real(bet))==0];
[bet,alpha]=solve(eqns,[bet,alpha]);
B=[-kappa,-(delta+2*g0*real(bet)),-2*g0*imag(alpha),0;...
delta+2*g0*real(bet),-kappa,2*g0*real(alpha),0;...
0,0,-gamma,omega_m;...
2*g0*real(alpha),2*g0*imag(alpha),-omega_m,-gamma;];
C=[0,-J,0,0;...
J,0,0,0;...
0,0,0,-K;...
0,0,K,0;];
H=[0,0,0,0;...
0,0,0,0;...
0,0,0,0;...
0,0,0,0;];
v