پایتون‌کار

from mpmath import *
# mpmath Option
mp.dps = 8; mp.pretty = True

#------ Input ----------------
def f(x,y):
    return x**2 + 2*x + y**2 - 8
def g(x,y):
    return x**2 - 2*x*y - 4
ep = 10**(-6)
x0 = 2
y0 = 4

#------ calc -----------------

def fx(x,y):
    return diff(f, (x, y), (1,0))
def fy(x,y):
    return diff(f, (x, y), (0,1))

def gx(x,y):
    return diff(g, (x, y), (1,0))
def gy(x,y):
    return diff(g, (x, y), (0,1))

def h0(x,y):
    return( g(x,y) * fy(x,y) - f(x,y) * gy(x,y) ) / \
          ( gy(x,y) * fx(x,y) - fy(x,y) * gx(x,y) )
def k0(x,y):
    return( g(x,y) * fx(x,y) - f(x,y) * gx(x,y) ) / \
          (-gy(x,y) * fx(x,y) + fy(x,y) * gx(x,y) )
    
xn=[x0]+[0 for i in range(20)]
yn=[y0]+[0 for i in range(20)]

xn[1] = xn[0] + h0(xn[0],yn[0])
yn[1] = yn[0] + k0(xn[0],yn[0])

n=1;
while ( abs( yn[n] - yn[n-1] ) > ep ):
    n+=1
    xn[n] = xn[n-1] + h0(xn[n-1],yn[n-1])
    yn[n] = yn[n-1] + k0(xn[n-1],yn[n-1])
    
print("n\tx[n]\t\t\ty[n]")
[print("{:}\t{:<{width}}\t{:<{width}}".format(i,xn[i],yn[i],width=16)) \
 for i in range(0,n+1)]

nst = len(str(round(1/mpf(ep))))-1

#------ Print -----------------
print('\nResult (x, y) = ({:}, {:})'.format(nstr(xn[n],nst),nstr(yn[n],nst)))

n	x[n]			y[n]
0	2               	4               
1	-6.0            	8.0             
2	-2.6829268      	5.0731707       
3	-1.221335       	3.6293203       
4	-0.70455505     	3.0793286       
5	-0.61300634     	2.978065        
6	-0.60986169     	2.9745276       
7	-0.60985793     	2.9745234       
8	-0.60985793     	2.9745234       

Result (x, y) = (-0.609858, 2.97452)

همواره دقت محاسبات باید بیشتر از دقت جواب نهایی باشد.