% compute and plot bifurcation diagrams

clear
global a c
a=0.2;
clf;
for c0=-0.15:0.005:0.15
    c=c0;
    for u0=-0.1:0.1:1.1
        ustar=fzero('fcubic',u0);
        plot(c,ustar,'.');
        hold on
    end;
end;
xlabel('c','fontsize',16)
ylabel('u*','fontsize',16)

pause;

u=[-0.1:0.01:1.1];
cstar2=u.*(u-a).*(u-1);
plot(cstar2,u,'g')

pause;

figure(2)
clf
c=0.0;
clf;
for a0=-0.15:0.01:0.35
    a=a0;
    for u0=-0.1:0.1:1.1
        ustar=fzero('fcubic',u0);
        plot(a,ustar,'.');
        hold on
    end;
end;
axis([-0.15,0.35,-0.2,1.2]);
xlabel('a','fontsize',16)
ylabel('u*','fontsize',16)

pause;

astar1=u;
plot(astar1,u,'g');
plot([-0.15,0.35],[0,0],'g');
plot([-0.15,0.15],[1,1],'g');
axis([-0.15,0.35,-0.2,1.2]);
