function res = mynewton(fun,fun_der,x0,err,n_steps)
% res = mynewton(fun,fun_der,x0,err,n_steps) tries to find a solution of
% the equation fun(x) = 0 using Newton's method with starting value x0. The
% code is intended for functions in one variable, although it is possible
% that it works also for functions in n variables (this is not tested!!!).
% The input fun_der should be the exact derivative of fun, else the
% convergence order is most probably not quadratic.
%
% If err > 0 is given, the computations stop as soon as the updates become
% (in norm) smaller than err; else they stop when they are smaller than
% machine accuracy.
% If n_steps is given, the code performs at most n_steps iterations; else
% n_steps is set to be 50.


% set the default number of iterations to 50
if(nargin < 5)
    n_steps = 50;
end

% set the default accuracy to 10^(-12)
if(nargin < 4)
    err = eps;
end

% Initialize the result
res = x0;

% Main iteration
for k=1:n_steps
    % compute the update step
    step = -feval(fun_der,res)\feval(fun,res);
    
    % check whether the update step is ok
    if(isnan(step) || isinf(step))
        error('Iteration stopped because of division by zero.');
    end
    
    % update of the iterate
    res = res + step;

    % stop the iteration, if the updates become too small
    if(norm(step) < err)
        disp(['Iteration converged after ',num2str(k),' steps.']);
        break;
    end
end