% NEWTEXAMPLE
% A function to show how sensitive Newton's Method is to the initial guess.
% 
% Syntax: newtexample(poly,[a,b],N,tol,maxits)
%
% Inputs: poly -- a vector holding the coefficients of a polynomial, from
%                 highest-order term to lowest (constant).  Use zero as a
%                 placeholder for any missing terms.  E.g. [2 1 0 -5]
%                 represents the polynomial 2x^3 + x^2 - 5.
%         [a,b] -- a vector giving the range of x values to use as starting
%                  values (x0) for Newton's Method.
%         N -- [optional] number of x values to use within the range [a,b]
%              (including a and b themselves).  Default = 101.
%         tol -- [optional] tolerance to use for Newton's Method.
%                Iterations will stop once |f(x)| < tol OR |dx| < tol,
%                where dx = change in x from one iteration to the next.
%                Default = 10^(-8).
%         maxits -- [optional] maximum number of iterations to perform
%                   before deciding that there is no convergence.
%                   Default = 50
%
% Outputs: no numeric output; makes three graphs:
%          1) graph of function over [a,b]
%          2) graph of number of iterations required to converge for each
%             starting x value (x0)
%          3) graph of final x value after convergence for each starting x
%             value (x0); values for which there was no convergence are not
%             shown
%
% Examples:
% >> newtexample([2,-3,-2],[-2,4])
%    p(x) = 2x^2 - 3x - 2 over domain [-2,4], with default options
%
% >> newtexample([2,-3,-2],[-2,4],501)
%    same as above but with more x points
%
% >> newtexample([2,-3,-2],[-2,4],501,1e-10,100)
%    same as above but with more points, a stopping tolerance of 10^(-10),
%    and a maximum iteration limit of 100
%
% >> newtexample([-1,3,0,-2],[-2,4],501)
%    p(x) = -x^3 + 3x^2 - 2 over [-2,4] with 501 points

% function definition
function newtexample(c,xr,N,tol,maxit)

% Set defaults for optional arguments
if (nargin<5),
    maxit = 50;
    if (nargin<4),
        tol = 1e-8;
        if (nargin<3),
            N = 101;
        end
    end
end

% find degree of polynomial (actually degree - 1)
ndeg = length(c)-1;
% calculate coefficients of derivative
dc = (ndeg:-1:1).*c(1:ndeg);

% set up array of x values and set aside space for answers
x = linspace(xr(1),xr(2),N);
its = zeros(size(x));
rts = its;

% loop through x values
for j=1:N
    thisx = x(j);  % get starting x value
    k = 0;  % initialize counter
    fval = polyval(c,thisx);  % evaluate function
    dx = 1+tol;  % no dx yet so define it to be greater than tol
    % loop until convergence -- ie while NOT converged and NOT over max
    % number of iterations
    while ((abs(fval)>tol) && (abs(dx)>tol) && (k < maxit)),
        dx = fval/(polyval(dc,thisx));  % calculate change
        thisx = thisx - dx;  % update x
        fval = polyval(c,thisx);  % re-evaluate function
        k = k+1;  % another iteration done (update counter)
    end
    % why did we quit?  convergence or b/c max iterations reached?
    if (k < maxit)
        % convergence, so save which x value we ended up at
        rts(j) = thisx;
    else
        % didn't converge, so x value = garbage (set = infinity)
        rts(j) = Inf;
    end
    % save the number of iterations made
    its(j) = k;
end

% show results
figure(1)
subplot(311)  % subdivide plot region into three
plot(x,polyval(c,x),x,0*x,'--')  % plot the function itself
% create string for graph title
tstr = ['graph of function: '];
for j=1:ndeg
    tstr = [tstr,num2str(c(j)),'x^',num2str(ndeg-j+1),' + '];
end
tstr = [tstr,num2str(c(ndeg+1))];
% label the first graph
xlabel('x'), ylabel('f(x)'), title(tstr)

subplot(312)  % next plot
plot(x,its,'.')  % plot number of iterations
% labels
tstr = 'Number of iterations for Newton''s method to converge.  Max = ';
tstr = [tstr,num2str(maxit)];
xlabel('x_0'), ylabel('iterations'), title(tstr)

subplot(313)  % next plot
plot(x,rts,'.')  % plot which roots were found
% labels
tstr = 'Root to which Newton''s method converges';
xlabel('x_0'), ylabel('root found'), title(tstr)