clear all


epsilon =-2;
h = +.4;

GAMMA = linspace(-20,15,201);

for ig=1:length(GAMMA)

gamma = GAMMA(ig);

C = [-1 0 -gamma 0 epsilon -h];

R = roots(C);

clear Xe
iter=1;
for ir=1:5
    if abs(imag(R(ir)))<1e-12
        Xe(iter) = R(ir);
        iter=iter+1;
    end
end



figure(1)
plot(gamma*ones(length(Xe),1),Xe,'.b')
hold on;
drawnow;
title('bifurcation diagram','fontsize',20)
xlabel('gamma','fontsize',20)
ylabel('x_e','fontsize',20)
end