newton法を用いた複素関数の解の導出

Question

ニュートン法を用いた複素関数の解の導出(Fortran)

本サイトの他質問者の投稿を参考にして以下のプログラムを組みました。
しかし、コンパイルは通るのですが、実行ができず、画面が止まったままになります。
初期値をいじると、x0のみ出力されることはありますが、
想像している解では無く困っております。
誰か分かる方教えてください。

発生している問題・エラーメッセージ

x0=( 0.500 0.500), ans=( 0.500 0.500)

このまま止まっています。

該当のソースコード

Implicit none
complex(kind(0d0)),parameter:: &

  x0=( 0.5, 0.5), &
  x1=( 1.0, 1.0), &
  x2=( 0.5, 0.0)
  write(*,'("x0=(",2f7.3,"), ans=(",2f7.3,")")') x0,newton(x0)
  write(*,'("x1=(",2f7.3,"), ans=(",2f7.3,")")') x1,newton(x1)
  write(*,'("x2=(",2f7.3,"), ans=(",2f7.3,")")') x2,newton(x2)
  stop
contains

  function newton(x0)
    Integer, parameter :: N = 500
    Integer, parameter :: c = 50
    double precision, parameter :: r1 = 1.d-2
    double precision, parameter :: delta = 1.d-3

    complex(kind(0d0)) :: x0,xn, x,ff,df, newton
    ff(x)= cos(x*(N+1.d0))-cos(x*(N+1.d0-2.d0*c))-r1/delta*sin(x)*sin(x*(N+1.d0))
    df(x)= -(N+1.d0)*sin(x*(N+1))+(N+1.d0-2.d0*c)*sin(x)&
           -r1/delta*cos(x)*sin(x*(N+1.d0))+sin(x*(N+1.d0))*cos(x*(N+1))
    x=x0
    do while(.true.)
      xn=x - ff(x)/df(x)
      if(abs(xn - x) < 1d-10)exit
      x=xn
    end do
    newton=xn
  end function
end

試したこと

様々な初期値を試しましたが、実行することができませんでした。
解析的に整理すると、0<Re[x]<πの範囲に解があり、
それらは実数か純虚数であると思われます。

補足情報（FW/ツールのバージョンなど）

ここにより詳細な情報を記載してください。

Answer 1

ニュートン法だとすると df は ff の微分になっているはずですが、プログラムでは微妙にそうなっていません。それが問題かもしれません。

別に元のままでもいいですが、
do while(.true.) ... end do は do ... end do でいけます。

倍精度を double precision で宣言するのは廃止予定ですし real(kind(1.0d0)) のほうが統一が取れます。