Newton Raphson Method – Roots of Two Coupled Non-Linear Equations | FORTRAN 95

PROGRAM NR_SYSTEM
IMPLICIT NONE

REAL::ea=6,es=0.01,x0,x1,y0,y1,f,g,f_dx,f_dy,g_dx,g_dy
INTEGER::c=0

PRINT *, "============================================" 
PRINT *, "Program To Find Roots of Two Coupled Non-Linear Equations Using Newton Raphson Method-[BY - www.BottomScience.com]" 
PRINT *, "============================================"
PRINT *,'Initial approximation?'
read(*,*)x0,y0

DO WHILE(ea>es)
x1=x0-((f(x0,y0)*g_dy(x0,y0))-(g(x0,y0)*f_dy(x0)))/((f_dx(x0,y0)*g_dy(x0,y0))-(f_dy(x0)*g_dx(y0)))

ea=abs(((x1-x0)/x1)*100)

y1=y0-((g(x0,y0)*f_dx(x0,y0))-(f(x0,y0)*g_dx(y0)))/((f_dx(x0,y0)*g_dy(x0,y0))-(f_dy(x0)*g_dx(y0)))

ea=abs(((y1-y0)/y1)*100)

x0=x1
y0=y1
c=c+1

IF(c>50) EXIT
PRINT *,'Current roots',x1,y1
END DO

PRINT *,'FINAL ROOT IS',x1,y1
END PROGRAM

REAL function f(x1,y1)
REAL::x1,y1
f=(x1**2)+(x1*y1)-10

return
end function

REAL function g(x2,y2)
REAL::x2,y2

g=y2+(3*x2*y2**2)-57

return
end function

REAL function f_dx(x2,y2)
REAL::x2,y2
f_dx=(2*x2)+(y2)

return
end function

REAL function f_dy(x2)
REAL::x2

f_dy=(x2)

return
end function

REAL function g_dx(y2)
REAL::y2

g_dx=(3*y2**2)

return
end function

REAL function g_dy(x2,y2)
REAL::x2,y2

g_dy=1+(6*x2*y2)

return
end function