円と曲線の交点 #1

古流記法なFORTRANプログラム。
x-y直交座標系において原点を中心とする単位円と曲線のにおける交点を求める。
曲線のy座標は負の値をとらないので、単位円（の半分）はと書ける。
求めるべきはとなるxである。
普通にニュートン・ラフソン法を使う。
微分係数は導関数を与えるのではなく2階中心差分で求める。

      PROGRAM INTERS
      IMPLICIT NONE
      DOUBLE PRECISION Y,YD,Y1,Y2,X,XN,E,H,X0
      PARAMETER(E=1.0D-12,H=1.0D-8,X0=0.9)
      X=X0
      DO
      XN=X-Y(X)/YD(X,H)
      IF(ABS(Y(XN)).LT.E)THEN
      WRITE(UNIT=*,FMT=*)'X=',XN
      WRITE(UNIT=*,FMT=*)'X*TAN(X)=',Y1(XN)
      WRITE(UNIT=*,FMT=*)'SQRT(1-X**2)=',Y2(XN)
      STOP
      END IF
      X=XN
      END DO
      END

      FUNCTION Y1(X)
      IMPLICIT NONE
      DOUBLE PRECISION Y1,X
      Y1=X*TAN(X)
      RETURN
      END

      FUNCTION Y2(X)
      IMPLICIT NONE
      DOUBLE PRECISION Y2,X
      Y2=SQRT(1-X**2)
      RETURN
      END

      FUNCTION Y(X)
      IMPLICIT NONE
      DOUBLE PRECISION Y,Y1,Y2,X
      Y=Y1(X)-Y2(X)
      RETURN
      END

      FUNCTION YD(X,H)
      IMPLICIT NONE
      DOUBLE PRECISION YD,Y,X,H
      YD=(Y(X+H)-2*Y(X)+Y(X-H))/H**2
      RETURN
      END

C     INTERSECTION POINT = (0.739085133,0.673612029)