function itstat = compndam ( dflag, x0, tol );
% CompNDam
%
% This m-file is part of the material of the lectures on
%   Numerical Mathematics
% at the Math Department of Halle University, Germany.
%
% Author :      Martin Arnold, martin.arnold@mathematik.uni-halle.de
% Version of :  Jan 16, 2006
%
% Computation with Damped Newton method
%
% Parameters:
%   dflag  (input)  : control flag "damping strategy"
%                       =0 .. none
%                       =1 .. brute force
%                       =2 .. Armijo like
%   x0     (input)  : initial guess
%   tol    (input)  : error tolerance (absolute and relative)
%   itstat (output) : iteration statistics ( NaN indicates values that are not applicable)
%                       [ k l xk fk fpk lam xnew |fnew/fpk| |fk/fpk| ]
%
% Examples:
%   itstat = compndam ( 0, sqrt(3), 1.0e-10 );
%   itstat = compndam ( 1, sqrt(3), 1.0e-10 );
%   itstat = compndam ( 1, 10, 1.0e-10 );
%   itstat = compndam ( 2, 10, 1.0e-10 );

colcurmax = 20;

atol  = tol;
rtol  = tol;

itstat = zeros(1,9);
colcur = 0;

iit  = 0;
xcur = x0;

cflag = 0;
while ~cflag, 

    fcur  = f(xcur);
    fpcur = fp(xcur);
    pcur  = - fcur / fpcur; 

    colcur = colcur + 1;
    itstat(colcur,[1 3:5 9]) = [ iit xcur fcur fpcur pcur ];

    iit    = iit + 1;
    lit    = 0;
    lamcur = 1; 

    lflag = 0;
    while ~lflag,

        xnew = xcur + lamcur*pcur;
        pnew = abs(f(xnew)/fpcur);
        itstat(colcur,[2 6:8]) = [ lit lamcur xnew pnew ];

        if dflag==0,
            lflag = 1;         
        elseif dflag==1,
            if pnew<=abs(pcur),
                lflag = 1;
            end;
        elseif dflag==2,
            if pnew<=(1-lamcur/2)*abs(pcur),
                lflag = 1;
            end;
        end;

        if ~lflag,
            lit = lit + 1;
            lamcur = lamcur/2; 
            colcur = colcur + 1;
            itstat(colcur,[1 3:5 9]) = NaN;
        end;

    end; 

    if ( (abs(xcur-xnew)<atol+rtol*abs(xcur) ) | ( colcur>colcurmax ) ),
        cflag = 1;
    end;

    xcur = xnew;

end;

iit = iit + 1;
colcur = colcur + 1;
itstat(colcur,[1 3]) = [ iit xcur ];
itstat(colcur,4:9)   = NaN;

function func = f ( x );
func = atan ( x );

function func = fp ( x );
func = 1 / ( 1 + x^2 );