f1 = @(x,y) ((x^2)/9) + ((x*y)/3) + y^2 -1
f2 = @(x,y) (x^2) - 1 - y

df1x = @(x,y) (2*x)/9 + (y/3) 
df1y = @(x,y) -(x/3) + (2*y)

df2x = @(x,y) 2*x
df2y = @(x,y) -1

DF = @(x,y) [ df1x(x,y) df1y(x,y);df2x(x,y) df2y(x,y)]

z = [.1;.1];
a = [-2.5;2.5];
c = [-0.175;-0.975];
d = [0;-1];

for j=1:50
    b = [f1(z(1),z(2));f2(z(1),z(2))];
    A = DF(z(1),z(2));
    z = z - A\b;
end

for j=1:50
    b = [f1(a(1),a(2));f2(a(1),a(2))];
    A = DF(a(1),a(2));
    a = a - A\b;
end

for j=1:50
    b = [f1(c(1),c(2));f2(c(1),c(2))];
    A = DF(c(1),c(2));
    c = c - A\b;
end

for j=1:50
    b = [f1(d(1),d(2));f2(d(1),d(2))];
    A = DF(d(1),d(2));
    d = d - A\b;
end


x = -sqrt(12):0.01:sqrt(12);

y1 = sqrt(1-x.^2/12)-x/6;

y2 = -sqrt(1-x.^2/12)-x/6;

y3 = x.^2 - 1;

plot(x,y1,'b'); hold on; plot(x,y2,'b'); hold on; plot(x,y3,'r');

hold on; plot(z(1),z(2),'o'); plot(a(1),a(2),'o'); plot(c(1),c(2),'o'); plot(d(1),d(2),'o');
z
a
d
c

 


    