  function ShowQNewton()
% Illustrate the Quasi-Newton Idea for systems of nonlinear
% equations.

x0 = [1;1;1];
tol = .000000000001;
itMax = 20;
xstar = QNewton(@MyF,@MyJ,x0,tol,itMax);

  function xStar = QNewton(F,J,x0,tol,itMax)
% Finds a zero of the function F:Rn ---> Rn with starting vector x0.
% Uses the Quasi-Newton Method.
% Terminates as soon as norm(F(xStar))<=tol or after itMax steps, whichever
% comes first.
xc = x0;
Fc = MyF(xc);
% Use finite differences to get the first Jacobian. Its QR factorization will be
% updated thereafter...
Jc = MyJ(@MyF,xc,Fc);
[Qc,Rc] = qr(Jc);
its = 0;
while norm(Fc)>tol && its<itMax
    % Compute the step using the current QR facorization
    sc = -(Rc\(Qc'*Fc));
    xPlus = xc + sc;
    FPlus = F(xPlus);
    % Looking at the linear part of the Taylor expansion, the next Jacobian 
    % JPlus should satisfy  Fc = FPlus + JPlus*(xc - xPlus), i.e.,
    % JPlus*sc = FPlus - Fc = yc
    yc = FPlus - Fc;
    % The problem of minimizing || JPlus - Jc ||_{F} subject to the 
    % constraint that JPlus *sc = yc is solved by setting
    % JPlus = Jc + (1/(sc'*sc)*(yc - Jc*sc)*sc' 
    % We use UpdateQR to get the QR factorization of this matrix...
    [Qc,Rc] = UpdateQR(Qc,Rc,(yc - Qc*(Rc*sc))/(sc'*sc),sc);
    Fc = FPlus;
    xc = xPlus;
    its = its+1;
    fprintf('%3d  %10.3e   %20.15f  %20.15f  %20.15f\n',its,norm(Fc),xc)
end
xStar = xc;

  function Fc = MyF(xc)
% xc is a column 3-vector
% This function is zeroed by a point that is on three given elliposids
A1 =  [6  -7   2; -7  10  -3;  2  -3  1];
A2 =  [5  -1   1; -1  10  -1;  1  -1  1];
A3 =  [9   4  -2;  4   2  -1; -2  -1  1];
v =  [2;1;-1];
Fc = [(xc-v)'*A1*(xc-v)-1;...
      (xc-v)'*A2*(xc-v)-1;...
      (xc-v)'*A3*(xc-v)-1];
 
  function Jc = MyJ(F,xc,Fc)
% The Jacobian of the function F at xc. Fc = F(xc)  
n = length(xc);
Jc = zeros(n,n);
for k=1:n
   % The divided difference parameter...
   delta = sqrt(eps)*xc(k);
   ek = [zeros(k-1,1);1;zeros(n-k,1)];
   % Approximate the partial of F(x) with respect to x(k)
   Jc(:,k) = (F(xc+delta*ek)-Fc)/delta;
end