ニュートン法のプログラム(FORTRAN)

3時間もかかった。f=x^2-2の例。もっとスマートなやり方があれば教えてください。
コンパイルするときに、頼むエラーでないでくれ〜って思っていると心臓が早鐘を打つ。

program newton_method
implicit none
real x1,x2,acc,root
integer eflag

eflag=0
acc=10.0**(-6)
x1=1.0
x2=2.0
!read *,x1,x2
call newton(x1,x2,acc,root,eflag)
if (eflag==0) then
print *,'root=',root
else
print *,'error=',eflag
endif
end program newton_method

subroutine newton(x1,x2,acc,root,eflag)
implicit none
real x1,x2,acc
real root
integer eflag

integer imax
parameter(imax=100)
integer i
real x0,xn1,xn2,f1,fd1,f2,fd2

x0=(x1+x2)/2.0
call funcd(x0,f1,fd1)
if (f1==0.0) then
root=x0
return
endif
xn1=x0

!----------------ループ
do i=1,imax
call funcd(xn1,f2,fd2)
xn2=xn1-(f2/fd2)

if (abs(xn2-xn1)<=acc) then
root=xn2
return
endif

if (xn2>=x2) then
eflag=1
return
endif
if (xn2<=x1) then
eflag=1
return
endif

xn1=xn2
enddo
!----------------ループ終

eflag=2
root=xn2
end subroutine newton

subroutine funcd(x,f,fd)
implicit none
real x
real f,fd
f=x**2-2
fd=2.0*x
end subroutine funcd