% find local extrema by solving for critical points
% NewtonsMethodSystems(1,0,1e-6,1e-6)  (1.0954, -0.2739)
% NewtonsMethodSystems(-1,0,1e-6,1e-6) (-1.0954, 0.2739)
function soln=NewtonsMethodSystems(x0,y0,tolrho,tolalpha)

% F=inline('(x-y)./(2*x.^2+8*y.^2+3)');
% [x,y]=meshgrid(-2:.01:2,-2:.01:2);
% z=F(x,y);
% figure;surf(x,y,z)

newX=[x0;y0]-jacobian(x0,y0)\f(x0,y0);
normf0=norm(f(x0,y0));

count=1;
while norm(f(newX(1),newX(2)))>tolrho*normf0+tolalpha && count<20
    x0=newX(1);y0=newX(2);
    newX=[x0;y0]-jacobian(x0,y0)\f(x0,y0);
    count=count+1;
end
soln=newX;

end

function z=f(x,y) % need to solve dF/dx=0 and dF/dy=0

z=[-2*x.^2+8*y.^2+4.*x.*y+3; -2*x.^2+8*y.^2-16*x.*y-3];

end

function J=jacobian(x,y)

J=[-4*x+4*y 16*y+4*x; -4*x-16*y 16*y-16*x];

end