% A script to make web diagrams
%
% The system is of the form x(k+1) = f(x(k))
%
% The function "f" is defined in a file called "rhs.m"
%

disp('Web diagrams');

if ~(exist('FNAME'))
	FNAME = 'rhs';
end;


if (length(FNAME)==0), FNAME='rhs'; end;
entry = input(['Enter name of the function (', FNAME ,')-> '],'s');
if (length(entry) == 0), entry = FNAME; end;
FNAME = entry;





x0 = input('Please enter the initial value -> ');
x0 = x0(1);
n = input('Please enter the number of iterations (default 10) -> ');
if (isempty(n)), n=10; end;


disp('Working.');

x = zeros(1,n+1);
x(1) = x0;
for k=1:n
	x(k+1) = feval(FNAME,x(k));
end;


disp('Computing graph.');


xmin = min(x);
xmax = max(x);
dx = xmax - xmin;
t = xmin-.1*dx : dx/100 : xmax+.1*dx;
z = feval(FNAME,t);
plot(t,t,'b',t,z,'r');

xmin = min([t,z,xmin]);
xmax = max([t,z,xmax]);
axis([xmin xmax xmin xmax]);
axis('square');

hold on;

xx = [x;x];
xx = xx(:);
xx = xx';


% yy = xx(1:length(xx)-2);
% yy = [0,yy];
% xx = xx(1:length(xx)-1);

yy = xx(2:length(xx));
yy(1) = 0;
xx = xx(1:length(xx)-1);


plot(xx,yy);
grid;
hold off;