%                   Function Mfile   Newton4 
%          Newton's method for dimensions 1,2,3 and 4.
%   The advantage of this code is that one can use any names for
%   variable, e.g. A, B,C instead of x y z.
%   The input arguments are character strings for the independent
%   variables and functions. 
%   Example:   If we want to solve system 
%              s-2+ lambda^2 = 0
%              cos(s) -lambda = 0,
%   we first write the character strings:
%        >> f1 = 's-2+lambda^2' ;
%        >> f2 = 'cos(s)-lambda';
%   It is important not to leave any spaces in the expressions for f1 and f2.
%   Then we make the call to newton4
%        >> newton4(f1,f2, 's', 'lambda')

%  User is then requested to provide a starting value. Code makes
%  five Newton iterations .

function out = newton4(varargin)
 
n = nargin/2;
if n == 1
   f1 = varargin{1};
   x1 = varargin{2} 
   F = inline(f1,x1);
   xx1 = sym(x1); ff1 = sym(f1);
   JJ = diff(ff1,xx1);
   J = inline(vectorize(JJ), x1);

   % begin the Newton loop
   p = input('enter the starting value   ')
   disp(' ')
   disp('   p         F(p)      error(Newton step)  ' )
   format compact
    [p, F(p)]

   for k = 1:5
      s = J(p)\F(p);
      p = p-s;
      [p,F(p), abs(s)]
   end

elseif n == 2
   f1 = varargin{1};
   f2 = varargin{2};
   x1 = varargin{3};
   x2 = varargin{4};
   Fchar = ['[', f1, ';',f2, ']'];
   F = inline(Fchar, x1,x2);
   xx1 = sym(x1); xx2 = sym(x2); ff1 = sym(f1); ff2 = sym(f2);
   JJ =jacobian([ff1, ff2], [xx1, xx2]);
   row1 = [char(JJ(1,1)),' ', char(JJ(1,2))];
   row2 = [char(JJ(2,1)),' ', char(JJ(2,2))];
   Jchar = ['[', row1, ';', row2, ']' ];
   J = inline(Jchar, x1,x2);

   % begin the Newton loop
   p = input('enter starting point as a col. vector   ')
   disp(' ')
   disp('  p      F(p)      error(Newton step)  ')
   format compact
   [ p, F(p(1),p(2))]

   for k = 1:5
      s = J(p(1),p(2))\F(p(1), p(2));
      p = p-s;
      [p, F(p(1), p(2)), [norm(s); 0] ]
   end  
   
elseif n == 3
   f1 = varargin{1}; f2 = varargin{2}; f3 = varargin{3};
   x1 = varargin{4}; x2 = varargin{5}; x3 = varargin{6};
   Fchar = ['[', f1, ';', f2, ';', f3,']'];
   F = inline(Fchar, x1,x2,x3)
   xx1 = sym(x1); xx2 = sym(x2); xx3 = sym(x3);
   ff1 = sym(f1); ff2 = sym(f2); ff3 = sym(f3);
   JJ = jacobian([ff1,ff2,ff3], [xx1,xx2, xx3]);
   row1 = [char(JJ(1,1)), ' ', char(JJ(1,2)), ' ', char(JJ(1,3))];
   row2 = [char(JJ(2,1)), ' ', char(JJ(2,2)), ' ', char(JJ(2,3))];
   row3 = [char(JJ(3,1)), ' ', char(JJ(3,2)), ' ', char(JJ(3,3))];
   Jchar = ['[',row1, ';', row2, ';', row3, ']'];
   J = inline(Jchar, x1,x2,x3);
   % begin the Newton loop
   p = input('enter the starting point as a col. vector   ')
   disp(' ')
   disp('      p           F(p)            error(Newton step) ')
   format compact
   [p, F(p(1),p(2),p(3)) ]
   for k = 1:5
      s = J(p(1),p(2),p(3))\F(p(1), p(2), p(3));
      p = p - s;
      [p, F(p(1), p(2), p(3)), [norm(s); 0; 0] ]
   end
elseif n == 4
   f1 = varargin{1}; f2 = varargin{2};
   f3 = varargin{3}; f4 = varargin{4};
   x1 = varargin{5}; x2 = varargin{6};
   x3 = varargin{7}; x4 = varargin{8};
   Fchar = ['[', f1, ';', f2, ';', f3, ';', f4,']'];
   F = inline(Fchar, x1,x2,x3,x4);
   xx1 = sym(x1); xx2 = sym(x2); xx3 = sym(x3); xx4 = sym(x4);
   ff1 = sym(f1); ff2 = sym(f2); ff3 = sym(f3); ff4 = sym(f4);
   JJ= jacobian([ff1,ff2,ff3, ff4], [xx1,xx2,xx3, xx4]);
   r1 = [char(JJ(1,1)), ' ', char(JJ(1,2)), ' ', char(JJ(1,3)), ' ', char(JJ(1,4))];
   r2 = [char(JJ(2,1)), ' ', char(JJ(2,2)), ' ', char(JJ(2,3)), ' ', char(JJ(2,4))];
   r3 = [char(JJ(3,1)), ' ', char(JJ(3,2)), ' ', char(JJ(3,3)), ' ', char(JJ(3,4))];
   r4 = [char(JJ(4,1)), ' ', char(JJ(4,2)), ' ', char(JJ(4,3)), ' ', char(JJ(4,4))];
   Jchar = ['[', r1, ';', r2, ';', r3, ';', r4, ']'];
   J = inline(Jchar, x1,x2,x3,x4);

   % begin the Newton loop
   p = input('enter the starting point as a col. vector   ')
   disp(' ')
   disp('   p           F(p)        error(Newton step)  ')
   format compact
   [p, F(p(1),p(2),p(3),p(4)) ]
   for k = 1:5
      s = J(p(1), p(2), p(3), p(4))\F(p(1), p(2), p(3), p(4));
      p = p-s;
      [p, F(p(1), p(2), p(3), p(4)), [norm(s); 0 ; 0 ; 0]]  
   end
 end
 format loose