function [res_broyden, tot_iter] = broyden(Func, Jacobian, p0, maxit, acc)

%Starting Evaluations
accur=0.5*10^(-acc);

% Initialisation
stop=false;
k=0;
X0 = p0';  % p0 is a row vector, so X0 is column vector 
res_broyden = X0';

% Calculate the Jacobian for X0
A=feval(Jacobian, X0);

% Step 1: Newton method for first approach
F0 = feval(Func, X0);
X1 = X0 - A\F0;         % left matrix devide, is the same as inv(A)*F0
k = k + 1;
F1 = feval(Func, X1);

res_broyden = [res_broyden ; X1'];
Diffx(k) = norm(X1-X0,2);  


% Step 2: Repeat until condition is satisfied
while (stop==false)
    
	% Estimate the Jacobian
    sk = X1 - X0;
    yk = F1 - F0;
    skT = sk';
    
    s1 = yk - A*sk;
    s2 = skT * sk;
    
    A = A + (s1 * skT)/s2;
    
    % Solve the linear system 
    Dx = -A\F1;
	
    % Calculate the next approach point
    X0 = X1;
    X1 = X0 + Dx;
	F1 = feval(Func, X1);
    k = k + 1;
    
    res_broyden = [res_broyden ; X1'];
    Diffx(k) = norm(X1-X0,2);
    
    % Check the criterion
	if (norm(F1,2)<=accur) || (norm(Dx,2)<=accur)
		stop=true;
    elseif (k > maxit)  
		stop=true;
    end
   
end

tot_iter = k;
%plot(Diffx)