function x = newton2(funk,x0,tol)
% newton2 - newton algorithm for the system of equations f(x)=0
% 
%   Syntax: 
%           x = newton2(f,x0,tol)
%   Arguments: 
%           f - string containing the name of a function file  
%           x0 - an nx1 matrix, the initial approximation
%           tol - a positive number, the tolerance 
%   Returns: 
%           x - an nx1 matrix, an approximate solution 
%   Description:
%           The program newton2 uses Newton's method to compute an 
%           approximate solution of the system of equations f(x)=0.  
%           The file f.m must contain the function y=f(x).  
%           The Jacobian matrix A=f'(x) is computed by the subfunction
%           jacobi.m below.  
%           Here x,x0,y are of type nx1 and A is of type nxn.
%           If the initial approximation x0 is sufficiently close 
%           to an exact root x_exact, then the program
%           computes an approximate solution x with the error
%           |x-x_exact| < tol.  
%   Examples: 
%          x = newton2('sin', 3, 1e-7) computes pi to 7 decimals
%          x = newton2('funk1', 3, 1e-7) computes the square 
%          root of 2
%   See also:
%          bisec.m, fixpoint.m, regula_falsi.m, newton2.m
%

%----------------------------------------------------------------------
%   Date created: 1999-12-05
%
%   Author:  Stig Larsson
%------------------------------+---------------------------------------+

x = x0;
h = tol + 1;

while norm(h)>tol
   A = jacobi(funk,x); 
   b = -feval(funk,x);
   h = A\b;
   x = x + h;
end

%--subfunction--------------------------------------------------------

function A=jacobi(funk,x)

% compute Jacobian matrix of funk at x

h = 1e-6; % stepsize in central difference approximation

A = zeros(size(x));

for i = 1:size(x,1)

   e = zeros(size(x));
   e(i) = 1;              % e is the i:th basis vector

   % compute i:th column of A 
   A(:,i) = ( feval(funk,x+h*e) - feval(funk,x-h*e) )/(2*h);

end


