title: How to solve system of (non-linear) equations using Jacobian and Newton's Method in Matlab
body: <p>I am trying to build a function that can solve a system of n-1 (non-linear) equations with n unknowns in Matlab making use of Newton's method.</p>
  
  <p>I did some research online and came to the following procedure:</p>
  
  <p>Use matlab's 'null' function applied to an array of first partial derivatives (Jacobian) to obtain an m-element vector, (call this vector v,) which is to be orthogonal to all the m-1 vectors of the partials. After proceeding along this vector a distance h0, use the 'null' function again, this time to find the m-1 dimensional subspace orthogonal to v from that advanced point. Projecting all the original coordinates to coordinates in this subspace you can perform an m-1-dimensional Newton-Raphson procedure which should converge very rapidly provided your step h0 is not too large. Then convert this back to the original coordinates.</p>
  
  <p>I have the following system of equations:</p>
  
  <pre><code>fun = @(x,y,z)[x.^3+y.^2+z.^2,x.^2-y.^3+sin(z)]
  </code></pre>
  
  <p>and it's Jacobian:</p>
  
  <pre><code>Jac = @(x,y,z)[3*x^2,    2*y,    2*z; 2*x, -3*y^2, cos(z)]
  </code></pre>
  
  <p>And an initial guess: 
      x0 = [1,1,1]</p>
  
  <p>In my function I convert the vector to numbers, so I can work with them:
      x = num2cell(x0)</p>
  
  <p>Then I want to do N iterations:</p>
  
  <pre><code>for numIterations = 1:N
  
  %calculate the function values at the current iteration
  F = fun(x{:})
  
  %calculate the jacobian matrix
  J = Jac(x{:})
  
  %find direction in which the curve extends:
  v = null(J)
  
  xnew = xold + h0 * v
  xnew = xnew'
  
  subspace = null(xnew)
  
  nextX = multinewton(f,df,xnew)
  </code></pre>
  
  <p>Now I am stuck in the part where I use Newton's method to calculate the next coordinates. How do I do this?</p>
  
  <p>The Newton method I want to use to calculate the next coordinates, is the following method:</p>
  
  <pre><code>function [zero,res,niter]=newton(f,df,x0,tol,nmax,varargin)
  %NEWTON Find function zeros.
  %   ZERO=NEWTON(FUN,DFUN,X0,TOL,NMAX) tries to find the zero ZERO of the 
  %   continuous and differentiable function FUN nearest to X0 using the    Newton 
  %   method. FUN and its derivative DFUN accept real scalar input x and returns 
  %   a real scalar value. If the search fails an error message is displayed.
  %   FUN and DFUN can also be inline objects.
  %
  %   [ZERO,RES,NITER]= NEWTON(FUN,...) returns the value of the residual in ZERO
  %   and the iteration number at which ZERO was computed.
  
  x = x0; 
  fx = feval(f,x,varargin{:}); 
  dfx = feval(df,x,varargin{:});
  niter = 0; 
  diff = tol+1;
  resultVector = [x];
  while diff &gt;= tol &amp;&amp; niter &lt;= nmax
  disp(x)
  niter = niter + 1;
  diff = - dfx\fx;
  x = x + diff;     
  diff = abs(diff);
  fx = feval(f,x,varargin{:});        
  dfx = feval(df,x,varargin{:});
  resultVector = [resultVector x];
  end
  if niter &gt; nmax
  fprintf(['newton stopped without converging to the desired tolerance',...
  'because the maximum number of iterations was reached\n']);
  end
  res = fx
  zero = x;
  disp('Number of iterations is:');
  disp(niter);
  
  plot(resultVector);
  return
  </code></pre>
  
  <p>I think I do not completely understand the part 'projecting all the original coordinates to coordinates in this subspace you can perform an m-1-dimensional Newton-Raphson procedure'. What is meant by this?
  What do I give as input to the newton method?</p>
  
  <p>Many thanks in advance!</p>