function xsol = newton2d(f1, f2, xvar, yvar, xstart, ystart)
%NEWTON2D 2-dimensional Newton's method
% The command xsol = newton2d(f1, f2, xvar, yvar, xstart, ystart)
% searches for a vector root of the (nonlinear) system of symbolic
% equations f1(xvar, yvar) = 0, f2(xvar, yvar) = 0, starting from the initial 
% guesses xvar = xstart, yvar = ystart.  Here f1 and f2 are symbolic
% expressions in xvar and yvar.  The output is a 2-dimensional vector
% of class "double".
%
% Example: xsol = newton2d(x^2 + y^2 - 4, x - y, x, y, 1, 1) 
% solves the simultaneous equations x^2 + y^2 = 4, x - y = 0
% starting at the point [1, 1] and gives the output 
% xsol =
%
%    [1.4142, 1.4142]
%
%
% written by Jonathan Rosenberg, 8/10/99
% revised 11/07/16

% Start by computing partial derivatives.
f11 = diff(f1, xvar);
f12 = diff(f1, yvar);
f21 = diff(f2, xvar);
f22 = diff(f2, yvar);
% Vector functions for evaluation
F1  = str2func(['@(',char(xvar),',',char(yvar),')', vectorize(char(f1))]);
F2  = str2func(['@(',char(xvar),',',char(yvar),')', vectorize(char(f2))]);
A11 = str2func(['@(',char(xvar),',',char(yvar),')', vectorize(char(f11))]);
A12 = str2func(['@(',char(xvar),',',char(yvar),')', vectorize(char(f12))]);
A21 = str2func(['@(',char(xvar),',',char(yvar),')', vectorize(char(f21))]);
A22 = str2func(['@(',char(xvar),',',char(yvar),')', vectorize(char(f22))]);
X  = [xstart; ystart];
for counter=1:20
A  = [A11(X(1), X(2)), A12(X(1), X(2)); A21(X(1), X(2)), A22(X(1), X(2))]; 
F = [F1(X(1), X(2)); F2(X(1), X(2))];
    if max(abs(F)) < eps, xsol=[X(1), X(2)];
% process has converged
    return; end
    if condest(A) > 10^10, error('matrix of partial derivatives singular or almost singular -- try again with different starting values'); end
    X  = X - A \ F;
    end
xsol=[X(1),X(2)];