function [sols,Niter] = mynewt(f,x0,derf)
% Find roots of given function f.
% [sols,Niter] = mynewt(f,x0,derf)
% Note: 
%  f  -- objective function
%  x0 -- vector of initial guesses
%  derf -- derivative of f, if not given, mynewt uses 
%  differential quotient written as a subfunction numdiff 
%                   
%  Values returned:
%  sols -- vector of roots found
%  Niter -- Nr. of iterations needed (Niter(k)>99 -> NaN)
%  NOTE: Don't use 'roots' as a name of a variable in command window,
%   since 'roots' is the name of a function (for polynomial roots).

%{
 1. function [sols,Niter] = mynewt(f,x0,derf)
    Write numdiff to allow 2 argument-calls
   

  Note! In this version Niter is scalar, goal is to get all
  zeros, not to analyze crazy starting values.

  2.  Modify code s.t. Niter is vector, N(k) is # of iter. 
   for starting component x0(k).
  
%}
Niter =0;
x =x0;
h = sqrt(eps);

while(max(abs(f(x)))>1e-8 && Niter<100)
  if nargin==2
      df = numdiff(f,x,h); % Subfunction numdiff below
  else 
      df = derf(x);
  end
    x = x - (f(x)./df);
    Niter = Niter+1;
end
sols = x;

% Subfunction:

  function y = numdiff(f,x,h)
  % Finite difference derivative of f
  y=(f(x+h)-f(x))/h;
  end

end
%% 

