% 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
%
% This version is faster than slobif, but you need to write
% the RHS function more carefully so that 
% x and p can hold vectors of values

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

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


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

% 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 -> ');

% we start the x's in the middle of the range
x = ones(1,pres)*(x0+x1)/2;

for k=1:warmup
	x = feval(RHSNAME,x,pvec);
end;

for k=1:running
	x = feval(RHSNAME,x,pvec);
	idx = ceil( (xres-1).*(x-x0)./(x1-x0) ) + 1;
	hits = find( (idx>=1) & (idx <= xres));
	
	for i = 1:length(hits)
		pict(hits(i), idx(hits(i)) ) = 2;
	end;
	
end;
	



	
image(pict');
colormap([0 0 0; 1 1 1]);
axis('image');
axis('xy');
axis('off');