function [x0, k] = newton_method(xk1, fcn, err, K)
 
delta = 1e2*eps;         % for differentiation purposes
xk    = 10*xk1/(1+err);  % error high enough so the while cycle is entered
k     = 1;               % first step
while (abs((xk1-xk)/xk)  >= err) && (k < K)
    xk  = xk1;
    fx  = feval(fcn, xk);
    dfx = (feval(fcn, xk+delta) - feval(fcn, xk-delta)) / (2*delta);
    xk1 = xk - fx/dfx;
    disp([k xk1 fx]);
    k = k + 1;
end
x0 = xk;