% This solves the system
%     x0^2+4*y0^2-4     =0
%     y0-x0^2+.4*x0+1.96=0
% It uses Newton's method for solving systems
% 
% Give the initial guesses
%
x0=input('x0= '); y0=input('y0= ');
%
% Initialize variables.  
% The variable 'diff' is used to look at the difference between
%         successive iterates.
% The vector variable 'save_rt' is used to save the iterates and the 
%         values of diff for later examination.
% The variable 'eps' is the desired error tolerance.
% The variable 'loop' counts how many iterates have been computed.
%
diff=1;
save_rt=zeros(30,3);
save_rt(1,1)=x0; save_rt(1,2)=y0;
eps=1.0e-12;
loop=1;
%
% Begin the main loop.
%
while diff > eps
% Calculate the function values.
  f0=x0^2+4*y0^2-4;
  g0=y0-x0^2+.4*x0+1.96;
% Calculate the derivative values for the Jacobian matrix
  fx=2*x0; fy=8*y0;
  gx=-2*x0+.4; gy=1;
% Set up the Jacobian matrix and solve for the Newton correction.
  jac=[fx fy; gx gy];
  rhs=-[f0 g0]';
  delta=jac\rhs;
% Calculate the infinity norm of the correction.
  diff=norm(delta,inf);
% Calculate the new iterate.
  x1=x0+delta(1); y1=y0+delta(2);
% Increase loop, save the new iterate, and re-initialize.
  loop=loop+1;
  save_rt(loop,1:3)=[x1,y1,diff];
  x0=x1; y0=y1;
end
%
% Display the iterates and the norms of the successive differences.
%
format long e;
save_rt=save_rt(1:loop,:)
