%Script file to show the dynamics of points under x^2+c for
%  differing values of c.  In each choice of c, we choose
%  several points to graph using cobweb analysis to get an
%  idea of what is happening.  Keep in mind that there is a
%  saddle node bifurcation at c = 1/4, then another bifurcation
%  occurs at c = -3/4, then another at c = -5/4
%
%  The program will halt after each "pause" command until any key
%   is pressed (for example, after line 28).

echo on

%
% Set 1:  These graphs are for -3/4 < c < 1/4 
%         (Choice: c=-1/4)
%
%The fixed points are pplus and pneg below:
c=-1/4;
f=inline('x.^2-1/4');
pplus=(1/2)*(1+sqrt(1-4*c));
pneg=(1/2)*(1-sqrt(1-4*c));

%Choose points in the intervals: >pplus, between -pplus and +pplus
% and less than -p+:
x=[-pplus-0.1, linspace(-pplus, pplus, 10), pplus+0.1,pplus+0.2];
s=length(x);
for j=1:s
    y=cobweb2(f,x(j),10,-2,2,-2,2);
    pause
end


%
% Set 2:  -5/4 < c < -3/4  (Choice:  c=-1)
%
% the pause commands halts execution until a button is pressed:
pause
clear, clf, clc

c=-1;
f=inline('x.^2-1');
pplus=(1/2)*(1+sqrt(1-4*c));
pneg=(1/2)*(1-sqrt(1-4*c));

%Choose points in the intervals: >pplus, between -pplus and +pplus
% and less than -p+:
x=[-pplus-0.1, linspace(-pplus, pplus, 10), pplus+0.1,pplus+0.2];
s=length(x);
for j=1:s
    figure(1)
    y=cobweb2(f,x(j),20,-2,2,-2,2);
    figure(2)
    plot(y)
    pause
    figure(1)
    clf
end

%
% Set 3:  Change c < -5/4 (Choice:  c= -1.35)
%
pause
clear, clf, clc

c=-1.35;
f=inline('x.^2-1.35');
pplus=(1/2)*(1+sqrt(1-4*c));
pneg=(1/2)*(1-sqrt(1-4*c));

%Choose points in the intervals: >pplus, between -pplus and +pplus
% and less than -p+:
x=[-pplus-0.1, linspace(-pplus, pplus, 10), pplus+0.1,pplus+0.2];
s=length(x);
for j=1:s
    figure(1)
    y=cobweb2(f,x(j),20,-2,2,-2,2);
    figure(2)
    plot(y)
    pause
    figure(1)
    clf
end

%
% Set 4:  Take c = -2
%
pause
clear, clf, clc

c=-2;
f=inline('x.^2-2');
pplus=(1/2)*(1+sqrt(1-4*c));
pneg=(1/2)*(1-sqrt(1-4*c));

%Choose points in the intervals: >pplus, between -pplus and +pplus
% and less than -p+:
x=[-pplus-0.1, linspace(-pplus, pplus, 10), pplus+0.1,pplus+0.2];
s=length(x);
for j=1:s
    figure(1)
    y=cobweb2(f,x(j),50,-2.5,2.5,-2.5,2.5);
    figure(2)
    plot(y)
    pause
    figure(1)
    clf
end

echo off