syms x y;

symui = x - (x.^3) / 3 - y;
symvi = x + 0.5 - 3*y;

symdudx = diff(symui,x);
symdvdx = diff(symvi,x);
symdudy = diff(symui,y);
symdvdy = diff(symvi,y);

ex = 100000;
ey = 100000;

error = 10^-6;

x0 = -4;
y0 = 0; 

while (error < abs(ex) || error < abs(ey))
    dudx = subs(symdudx,[x,y],[x0,y0]);
    dvdx = subs(symdvdx,[x,y],[x0,y0]);
    dudy = subs(symdudy,[x,y],[x0,y0]);
    dvdy = subs(symdvdy,[x,y],[x0,y0]);
    ui = subs(symui,[x,y],[x0,y0]);
    vi = subs(symvi,[x,y],[x0,y0]);
    
    xjared = x0 - (ui*dvdy - vi*dudy)/(dudx*dvdy - dudy*dvdx);
    yjared = y0 - (vi*dudx - ui*dvdx)/(dudx*dvdy - dudy*dvdx);
    
    ex = (xjared - x0)/(xjared);
    ey = (yjared - y0)/(yjared);
    
    x0 = xjared;
    y0 = yjared;    
end

disp(x0);
disp(y0);

