%  Scripts for MAT 331 October 4, 2017 
%  Corrected after class

% function to find zero by bisection

function [x,k]=bisect(f,a,b)
if sign(a)==sign(b)
    if sign(a)==0
        x=a;
        k=0;
        return;
    else
       return;
    end
    k=0;
    while abs(b-a) > eps*abs(b)
        x = (a + b)/2;
      if sign(f(x)) == sign(f(b))
         b = x;
      else
         a = x;
      end
      [a,b];
      k = k + 1;
    end   
 return;

%---------------------------------------------------------

% function to apply Newton's method to a function
% derivative must be provided  

function [y,k] = newton(f,fprime,x)
   k = 0;
   xprev=x+1;
   while abs(x - xprev) > eps*abs(x)
      xprev = x;
      x = x - f(x)/fprime(x);
      k = k + 1;
   end
   y=x;
return;

%---------------------------------------------

% function to apply secant method. Function and
% two starting points must be given

function [out,k] = secant(f,a,b)
   k=0;
   while abs(b-a) > eps*abs(b)
      c = a;
      a = b;
      b = b + (b - c)/(f(c)/f(b)-1);
      k = k + 1;
   end
   out=a;
return;

%---------------------------------------------------------

% function to apply  Newton's method  to a polynomial
% derivative is computed by this function

function [y,k] = poly_newton(p,x)
   p_prime=polyder(p);
   k = 0;
   xprev=x+1;
   while abs(x - xprev) > eps*abs(x)
      xprev = x;
      x = x - polyval(p,x)/polyval(p_prime,x);
      k = k + 1;
      if k>50  % stop if doesn't converge
          y=NaN;
          return;
      end % if
   end % while
   y=x;
return;

% try p=[1,0,0,1] cube roots of 1
% try p=[2,8,2,0,3,1]


%--------------------------------------------------------


% function to plot result of Newton's method. Colors plane 
% according to which root that initial value finds

% inputs
% p = polynomial in MATLAB format: coefficeints in decreasing order
% a = size box, picture will be of [-a,a]^2
% n = number of pixels along each edge of box 

function   plot_newton(p,a,n)

   all_roots=roots(p);
   delta=abs(2*a)/n;
   % save_root=zeros(n);
   for k=1:n
      for j=1:n
         z=(-abs(a)+2*abs(a)*j/n)+1i*(-abs(a)+2*abs(a)*k/n);
         [y,iter]=poly_newton(p,z);
          if isnan(y)
              save_root(k,j)=0;
          else
              [m,root_num]=min(abs(all_roots-y));
             save_root(k,j)=root_num;
          end
      end
      k
   end
   all_roots;
   save_root;
   figure;
   title('Domains of attraction for Newtons method', p)
   image(save_root,'CDataMapping','scaled')
   %axis image
   %colormap(flipud(jet(length(all_roots))))
   toc;
return;


%--------------------------------------------
% function to count and find zeros

function  [counter,savezero]=count_zeros(f,a,b,n)
   delta=(b-a)/n;
   t=[a:delta:b];
   counter=0;
   clear savezero;
   for k=1:n-1
      if f(t(k))*f(t(k+1)) <= 0
         counter=counter+1;
         savezero=fzero(f,[t(k),t(k+1)]);
      end
   end 
return;