rbm=0.00595;
rbp=0.01279;
%b=0.16;
%a=0.001;
a=M.rgrid(1);
b=M.rgrid(end);

Q=mean(M.Er(1,floor(length(M.zgrid)/2),end-30:end),3)/2*a;
%Q=mean(M.Er(1,100,end-30:end),3)/2*a;
omega_d=M.omepe^2/2*M.omece;
if( a~=0)
    omega_q=2*Q/(M.B0*(rbm^2));
else
    omega_q=0;    
end

Levaluated=10;
stability=zeros(Levaluated,1);
for l=1:Levaluated
    stability(l)=(-l*(1-omega_q/omega_d)*(1-(rbm/rbp)^2)*(1-(a/b)^2)+ ...
        2*(1+(a/b)^(2*l))-(1+(rbm/rbp)^(2*l))*((a/rbm)^(2*l)+(rbp/rbm)^(2*l)))^2 -...
        4*(rbm/rbp)^(2*l)*(1-(rbp/b)^(2*l))^2*(1-(a/rbm)^(2*l))^2;
    if(stability(l) <0)
        fprintf(2,"Diocotron mode l=%d is potentially unstable\n",l);
    end
end
figure
plot(1:Levaluated,stability)