function plot1fun(funname,a,b,tol,details)
% Help for plot1fun.m
% copyright 1998, D. Estep
%
% PURPOSE
%
% This function plots a user input function on [a,b] at points determined
% using an adaptive strategy that insures the change in the function
% values from one sample point to the next is smaller than a tolerance.
% The program starts with a nearly uniform mesh of sample points and then
% creates an adapted mesh by successive refinement.
%
% FUNCTION STATEMENT
%
% plot1fun(funname,a,b,tol,details)
%
% REQUIRED VARIABLES
%
%    funname:     contains the name of the function to be plotted in the
%                 form functionname(x) where functionname is the name of the
%                 function and x is the variable
%    a,b:         numbers defining the interval [a,b] on which the function is plotted
%
% OPTIONAL VARIABLES
%
%    tol:         the tolerance for the maximum possible change in the function
%                 values at neighboring sample points. 
%    details:     indicates whether the program should give information on
%                 the sample points.
%                 details = 0 means just show plot
%                 details = 1 means show plot, changes in function, and mesh sizes
%
% USAGE
%
% First create a m file named functionname.m containing the function to be 
% plotted with the form
%
%    function [y] = functionname(x)
%    y = give the formula in terms of input x;
%
%    For example, we create a file called quadratic.m containing the lines
%    function [y] = quadratic(x)
%    y = x^2;
%
% At the matlab prompt, define a variable containing the string functionname(x)
%
%    For the example above, we would define 
%    >> funname = 'quadratic(x)'
%
% Also define the values of a, b, tol, details, then call
% >> plot1fun(variable containing function name, a, b, tol, details)
%
%    In the example above, we would type
%    >> plot1fun(funname,-4,4,.01,0)
%    to get a plot on the interval [-4,4] using a mesh of sample points
%    in which the change in the function values is smaller than .01
%
% The last two arguments to plotfun can be left off, for example
% >> plot1fun(funname,-4,4) 
%
% You may also plot standard built-in matlab functions, for example
% >> funname = 'sin(x)'
% >> plot1fun(funname,0,10)
% plots sin(x) on [0,10].
% >> funname = 'sin(x)*cos(x)*x'
% >> plot1fun(funname,0,10)
% plots sin(x)*cos(x)*x on [0,10]
%
% MORE EXAMPLES
%
% The following function requires denser mesh for larger x
% run with a = -2, b = -a, tol = .5 to .001
%
% function [y] = quadratic(x)
% y=x^2;
%
%
% The following function has two regions requiring more points
% run with a = - 3, b = 3, tol = .01
%
% function [y] = expo(x)
% y=exp(-10*(x+.5)*(x+.5)) + exp(-10*(x-2)*(x-2));
%
% To see how the adaptive strategy can fail, plot the function above
% using a = 4, b = -4, tol = .01 after commenting out the line below that
% adds a random perturbation to the initial uniformly distributed 
% sample points
%
% The following function requires many, many points near 0:
% run with a = 0, b = 3, tol = .1 to .01
%
% function [y] = sino(x)
% y=sin(1/(x+.05));
%
% DETAILS
%
% The program uses two meshes: a current mesh and a new mesh.  Going from
% left to right, the difference in function values at successive nodes in
% the current mesh is checked. If this difference is smaller than tol, the
% new mesh points are defined to be the old mesh points.  If this difference
% is larger than tol, then the interval between the nodes is divided into
% uniform pieces and the new mesh is defined by using the nodes of the
% pieces.  The number of pieces is computed by assuming the function is 
% linear between the nodes of the current mesh.  After every interval in
% the current mesh is checked, the current mesh is reset to be the just
% computed new mesh.  If no refinements were needed, the iteration stops
% and if there was a refinement in some interval, there is another loop.

%initialize

% define the maximum number of adaptive iterations
  maxrefineiter=5;

% define the function to be plotted
f = inline(funname,'x');

% check the number of arguments and set defaults for the missing values
noargs = nargin;
if noargs == 3
   tol = .1 *(b-a);
   details = 0;
elseif noargs == 4 
   details = 0;
elseif  noargs < 3  &  noargs>5 
   disp('wrong number of inputs!')
   break
end

% define initial number of sample points and uniform sample interval
n = 11;
delta = (b-a)/(n-1);

% define the sample points. The first and last points align with a and b.
x(1) = a;
y(1) = f(x(1));
x(n) = b;
y(n)=f(x(n));
for i = 2: n-1
   x(i) = a+ delta*(i-1);
	% vary the internal sample points by 10% randomly to avoid
	% problems with sampling at symmetric points
   x(i) = x(i) + (rand-.5)*.2*delta;
	y(i)=f(x(i));
end

% set refinement level
refine = -1;

% loopflag=1 means do a refinement search, the first time we do one of course
loopflag=1;

while loopflag>0 & refine <= maxrefineiter
   
   % xnew, ynew are the new mesh points and function values
   % nnew is the number of mesh points in the new mesh
   
   % the first nodes in two meshes always agree
   xnew(1) = x(1);
   ynew(1) = y(1);
   nnew=1;
   
   % default: no more refinement iterations
   loopflag = -1;
   refine = refine+1;
   
   for i = 2:n
      
      % check if function difference is small enough
      if (abs(y(i)-y(i-1))< tol)
         
         % if so, then set new mesh point = old mesh point and go to next
         nnew=nnew+1;
         xnew(nnew) = x(i);
         ynew(nnew) = y(i);
         
      else
         
         % indicate that a refinement is needed in this iteration, so check
         % for one more iteration to be sure that the differences are small
         loopflag = 1;
         
         % assume that the function is linear between the two mesh points and
         % compute the necessary number of divisions to get the desired change
         nmid=ceil(abs(y(i)-y(i-1))/tol)+1;
         delta = abs(x(i)-x(i-1))/nmid;
         
         % now define the new mesh points up to the next to last one
         for k = 1:nmid-1;
            xnew(nnew+k) = x(i-1) + delta*k;
            ynew(nnew+k) = f(xnew(nnew+k));
         end
         
         % treat the last point to be sure that it lines up with old mesh
         xnew(nnew+nmid) = x(i);
         ynew(nnew+nmid) = y(i);
         nnew = nnew + nmid;
                  
      end   
   
	end
   
   % prepare for a new loop, set the old values to be the new values
   clear x;
   clear y;
   for i = 1: nnew
      x(i) = xnew(i);
      y(i) = ynew(i);
   end
   n=nnew;
   
end

if refine > maxrefineiter
   fprintf('The program exceeded %d refinement iterations\n',maxrefineiter)
   break
end


if(details==0)
   plot(x,y)
   graphof=strcat('Graph of ',funname);
   title(graphof)
   xlabel('x')
else
   figure
   subplot(1,3,1)
   hold on
	plot(x,y)
   plot(x,0,'r.')
   graphof=['Graph of ',funname];
   title(graphof)
   xlabel('x')
	hold off
	for i = 1: n - 1
   	dy(i) = abs(y(i+1)-y(i));
      dx(i) = abs(x(i+1)-x(i));
      xx(i) = (x(i+1)+x(i))/2;
	end
	subplot(1,3,2)
   plot(xx,dy)
   title('Changes in Values')
   xlabel('x')
	subplot(1,3,3)
   semilogy(xx,dx)
   title('Interval Sizes')
   xlabel('x')
	fprintf('plotting took %d refinement steps\n', refine)

end