function broyden( s1,s2,a, n, tor)
%Broyden Method for 2 nonlinear equations
%
%epsilon = 10^-6;
syms x y ;
fun1=sym(s1);
fun2=sym(s2);
f1a = subs(fun1,[x y],a);
f2a = subs(fun2,[x y],a);
fa=[f1a f2a] ;
fprintf('initial a= %10.6f  %10.6f   f(a)= %0.7g %0.7g\n\n',a ,fa)
%start broyden
A0=jacobian([fun1; fun2],[x y]);
fprintf(' n           X                          F(X)\n')

for j = 1: n
%solve Ja* h = -fa'
   Ja=subs(A0, [x y],a);
   h=Ja\(-fa') ;
   f1b=subs(fun1,[x y],a);
   f2b=subs(fun2,[x y],a);
   fb=[f1b f2b] ;
   a = a+h' ; %generate next point
  f1a = subs(fun1,[x y],a);
  f2a = subs(fun2,[x y],a);
  fa=[f1a f2a] ;
  z=fa-fb;
  A0=A0+(z'+fb')*h'/(h'*h);
  fprintf('%3d    %10.6f  %10.6f     %0.7g  %0.7g\n',j,a ,fa)  
  if norm(h) < tor
     fprintf('\nconvergence\n')
     return
  end
end  %for     
