function [x,y,f,cvg_hst,dfdx] = HPGfsolve ( fctn, x, c, delta, tolerance )
% [x,y,f,cvg_hst,dfdx] = HPGfsolve ( fctn, x, c, delta, tolerance )
% solve a nonlinear system of equations  f = fctn(x,consts)
% using Newton-Raphson method
% 
% INPUT
% =====
% fctn     function of the form [f,y] = fctn(x,c) 
%   x      initial guess values of indpendent variables                n x 1
%   c      other constants passed to to the function  fctn(x,c)        1 x m
%   d      the increment used for the numerical jacobian               1 x 1
%   t      the relative convergence tolerance                          1 x 1
%
% OUTPUT
% =====
%   x      column vector of the solution f(x) = 0                      n x 1
%   f      a column vector of function values at the solution ~ zero   n x 1
% cvg_hst  a record of x and norm(f) from iteration to iteration    2n+1 x 1
%  dfdx    the jacobian at the optimal solution
%
% see: http://en.wikipedia.org/wiki/Newton%27s_method

if nargin < 4, delta = 1e-3; end % perturbation factor for numerical Jacobian
if nargin < 5, tolerance = 1e-9; end   % convergence tolerance

verbose   = 0;                % show convergence information
MaxIter   = 10;              % maximum number of Newton-Raphson iterations
nfeval    = 1;                % counter for the number of function evaluations

[nr,nc] = size(x);
n = nr*nc;                    % number of independent variables

[f,y] = feval( fctn, x, c );  % evaluate the function at variable val's "x"
norm_f = norm(f);             % the square root of the sum of the squares

m = length(f);

dfdx = zeros(m,n);                       % initialize Jacobian
cvg_hst = NaN(n+5,MaxIter);

for iter = 1:MaxIter                     % start Newton-Raphson iterations

   for j=1:n                             % compute the Jacobian numerically
 
     xp = x;
%    pxj = max(delta, abs(xp(j))*delta);    % perturbation in x(j)
     pxj = delta;
     xp(j) = xp(j) + pxj;                   % perturbed x(j)

     fp = feval( fctn, xp, c );             % function at perturbed x(j)

     dfdx(:,j) = (fp - f) / pxj;            % df/dx(j) . col. vector
 
   end                                   % end Jacobian computation

  dx = -dfdx \ f;                        % Newton-Raphson update direction
  dx = reshape(dx,nr,nc);

  [f_new,y_new] = feval( fctn, x+dx, c );
  norm_f_new = norm(f_new);

  nfeval = nfeval + n + 1;

  while norm_f_new > norm_f              % check for improvement
    dx = dx / 2.0;                       % reduce dx until solution improves
    [f_new,y_new] = feval( fctn, x+dx, c );
    norm_f_new = norm(f_new);
    nfeval = nfeval + 1;
  end
  x = x+dx;                              % update the independent variables
  cvg_hst(:,iter) = [ x; norm_f; -1; nfeval; norm(dx); norm(f-f_new) ]; % history

  f      = f_new;                        % update the function values
  norm_f = norm_f_new;                   % update the function norm
  y      = y_new;                        % update ohter outputs

  if verbose && iter > 1
    disp([ iter f_new'   x' norm_f  tolerance ]);
    %disp([ eig(dfdx) , f , x ]);
  end

  if  norm_f  < tolerance                % check for convergence
    break
  end

end                                      % end Newton-Raphson iterations
cvg_hst = cvg_hst(:,1:iter);
if (iter == MaxIter)
  fprintf(' !! HPGfsolve failed to converge to tolerance %e in MatIter = %d iterations!! \n', tolerance, MaxIter);
  norm_f
end

% HPGfsolve - HP Gavin - 19 Apr 2010, 4 Jan 2011, 17 Apr 2017, 
