% BROYDEN: Denne lo'ser f(x) med Broydens metode
% Kall: root=broyden_secant(func,x0,B)

function [root,B,iter1]=broyden_secant(func,x0,B)

global n x_pred x_old teller

nn=max(size(x0));
delta=0.0000001;
x=x0;
size(x);
f=feval(func,x);
size(f);
error=norm(f);
diff_f=norm(f);

if ~exist('B')
  B=[];
end

%Beregn Jacobi-matrise:
if isempty(B),
  teller=1;
  for i=1:nn
    p=0*x; p(i)=delta;
    J(:,i)=(feval(func,x+p)-feval(func,x-p))/(2*delta);
  end
  Js=sparse(J);
  %B=inv(J);
  B=inv(Js);
end %if

iter=0;
while error>1.e-6     %Main loop

%if rem(teller,20)==0
%  for i=1:nn
%    p=0*x; p(i)=delta;
%    J(:,i)=(feval(func,x+p)-feval(func,x-p))/(2*delta);
%  end
%
%  B=inv(J);
%end %if
  
  
  
  iter=iter+1;
  dx=-B*f;
  xnew=x+dx;
  if isnan(xnew)
      disp('dsfjk')
  end
  
  fnew=feval(func,xnew);

  count=0;
  while norm(fnew)>error  %Backtrack if norm(f) not decreasing
    count=count+1;
    if count>100; break;end;
    dx=dx/2;
    if count==4    %If not success on 4th backtrack, reevaluate Jacobian
      for i=1:nn
        p=0*x; p(i)=delta;
        J(:,i)=(feval(func,x+p)-feval(func,x-p))/(2*delta)';
      end
      B=inv(J);
      dx=-B*f;
    end %if
    xnew=x+dx;
   
    fnew=feval(func,xnew);
  end %while



df=fnew-f;
  temp=B*df;
  B=B+(dx-temp)*((dx'*B)/(dx'*temp));
  
  x=xnew;
  f=fnew;
  error=norm(f);

  if iter>100
    'IKKE KONVERGENS I BROYDEN'
     break
  end
end %while


root=x;
iter1=iter; 





