function [f,x] = NRNDN(x0,tol,iprint)
%
% function [f,x] = NRNDN(x0,tol,iprint)
% This function performs multivariate Newton-Raphson
% with numerical approximations to the partial derivatives
%
% x0 is the vector of initial guesses
% tol is the absolute tolerance on each function
% iprint = 1 if you want iteration information
%
% the functions must be located in syseqninput.m
%
% author: David Keffer, University of Tennessee, Knoxville
% date: May, 1999
%
%
% maximum iterations
maxit = 100;
n = max(size(x0));
dxcon = zeros(n,1);
dxcon(:) = 0.01;
Residual = zeros(n,1);
dx = zeros(n,1);
Jacobian = zeros(n,n);
xp = zeros(4);
fp = zeros(n,n,4);
%
% initialization
%
x = x0;
err = 100.0;
iter = 0;
%
% convergence loop
%
while ( err > tol )
for j = 1:1:n
dx(j) = min(dxcon(j)*x(j),dxcon(j));
i12dx(j) = 1.0/(12.0*dx(j));
end
%
% evaluate function
%
Residual = syseqninput(x);
Residual = Residual';
%
% set up a grid for evaluation of numerical partial derivatives
% evaluate function at grid points
%
for j = 1:1:n
   for i = 1:1:n
      xeval(i) = x(i);
end
xp(1) = x(j) - 2*dx(j);
xp(2) = x(j) - dx(j);
xp(3) = x(j) + dx(j);
xp(4) = x(j) + 2*dx(j);
for k = 1:1:4
xeval(j) = xp(k);
fp(:,j,k) = syseqninput(xeval);
end
end
%
% calculate numerical derivative
%
for i = 1:1:n
for j = 1:1:n
Jacobian(i,j) = i12dx(j)*( -fp(i,j,4) + 8*fp(i,j,3) - 8*fp(i,j,2) + fp(i,j,1) );
end
end
%
% evaluate new estimate of solution
%
xold = x;
invJ = inv(Jacobian);
deltax = -invJ*Residual;
for j = 1:1:n
x(j) = xold(j) + deltax(j);
end
%
% check for convergence
%
iter = iter +1;
err = sqrt( sum(deltax.^2) /n ) ;
if (iprint == 1)
fprintf (1,'iter = %4i, err = %9.2e \n ', iter, err);
end
if ( iter > maxit)
Residual
error ('maximum number of iterations exceeded');
end
end
f = sqrt(sum(Residual.*Residual)/n);