% A script to prepare a bifurcation diagram
% The function is of the form y = f(x,p)
% where x is the argument of the function and p is the parameter

global RHSNAME
disp('Bifuctation digram');
setrhs;

% the resolution of the image
% 90 x 90 should work in the student version
pres = 50;
xres = 40;


% warmup is the number of iterations to perform 
% before plotting. running is the number of points to 
% really compute.
warmup = 50;
running = 75;

% pict stores the picture
pict = ones(pres,xres);

% get parameter range:
p0 = input('Please enter smallest parameter value -> ');
p1 = input('Please enter largest parameter value -> ');

pvec = [p0 : (p1-p0)/(pres-1) : p1];

% get variable range:
x0 = input('Please enter smallest x value we should see -> ');
x1 = input('Please enter largest x value we should see -> ');
xmid = (x0+x1)/2;


for pidx = 1:pres		% for all values of the parameter...
	if (rem(pidx,10)==0) 
		disp(pidx)
	end;
	p = pvec(pidx);
	x = xmid;
	xvec = ones(1,xres);
	%  run warmup iterations 
	for k=1:warmup, x = feval(RHSNAME,x,p); end;
	
	% produce the values
	for k=1:running
		x = feval(RHSNAME,x,p);
		idx = ceil( (xres-1)*(x-x0)/(x1-x0) ) + 1;
		if (1 <= idx) & (idx <= xres)
			xvec(idx)= 2;
		end;
	end;
	
	% save the points in pict
	pict(pidx,:) = xvec;
end;
	
image(pict');
colormap([0 0 0; 1 1 1]);
axis('image');
axis('xy');
axis('off');