Newton Raphson Solver with adaptive Step Size

NewtonRaphson solves a system of nonlinear equations via Newtons method

% NewtonRaphson solves equations of the form:
%
% F(X) = 0 where F and X may be scalars or vectors
%
% NewtonRaphson implements the damped newton method with adaptive step
% size. Theory and discussion can be found here:
% http://forum.vis.ethz.ch/showthread.php?13434-Damped-Newton-method-Schrittweitensteuerung
%
% The Jacobian is calculated numerically by a simple forward differential
% method. Find more here:
% http://en.wikipedia.org/wiki/Numerical_differentiation
%
% The error estimation ist made by root-mean-square
% http://en.wikipedia.org/wiki/Root_mean_square
%
%
% x = NewtonRaphson(FUN,X0) starts at the initial guess X0 and tries to
% solve the equations in FUN. FUN is a function handle and has to accept
% input x and return a vector of equation values F evaluated at x.
% Default values for solver and display setting.
%
% x = NewtonRaphson(FUN,X0,lambda) starts at the initial guess X0 and tries to
% solve the equations in FUN with user supplied initial relaxation factor. Might
% be useful to increase solution speed.
% default value: lambda = 0.1
%
% x = NewtonRaphson(FUN,X0,lambda,maxIter) ...
% User supplied maximum number of iterations
% default value: maxIter = 100
%
% x = NewtonRaphson(FUN,X0,lambda,maxIter,Display) ...
% Display options:
% 'on' - full output during solution process
% 'off' - hide output
% default value: maxIter = 'off'
%
% OUTPUT Arguments:
% x - Solution
% F - Value at x
% Jac - Jacobian at x
% Exitflat: exit conditions of damped newton
% 1: all ok. Solution found
% -1: no solution found
% -2: FUN is not a function handle
%
% To enter Demo Mode start without any arguments
% NewtonRaphson();
%
% Examples:
%
% Find zero of function atan(x)
%
% create function handle (here as an anonymous function)
% F = @(x)atan(x);
% x0 = 5; % start value /initial guess
%
% x = NewtonRaphson(F,x0,1);
%
% Find solution of following system of equations (same as in fsolve help)
% F = @(x)[2*x(1) - x(2) - exp(-x(1));
% -x(1) + 2*x(2) - exp(-x(2))];
% x0 = [1;2];
%
% x = NewtonRaphson(F,x0,0.1,100,'on');