方程如下：

算法用的是徐士良的《Fortran常用算法程序集》中的二分法。

SUBROUTINE DDHRT(A,B,H,EPS,X,N,M,F)

DIMENSION X(N)
M=0
Z=A
Y=F(Z)

Y=F(Z)
GOTO 10
END IF
Z1=Z+H
Y1=F(Z1)
IF (ABS(Y1).LT.EPS) THEN
M=M+1
X(M)=Z1
Z=Z1+H/2.0
Y=F(Z)
GOTO 10
END IF
IF (Y*Y1.GT.0.0) THEN
Y=Y1
Z=Z1
GOTO 10
END IF
20 IF (ABS(Z1-Z).LT.EPS) THEN
M=M+1
X(M)=(Z1+Z)/2.0
Z=Z1+H/2.0
Y=F(Z)
GOTO 10
END IF
Z0=(Z1+Z)/2.0
Y0=F(Z0)
IF (ABS(Y0).LT.EPS) THEN
M=M+1
X(M)=Z0
Z=Z0+H/2.0
Y=F(Z)
GOTO 10
END IF
IF (Y*Y0.LT.0.0) THEN
Z1=Z0
Y1=Y0
ELSE
Z=Z0
Y=Y0
END IF
GOTO 20
END
c=========================

DIMENSION X(3)
EXTERNAL F
CALL DDHRT(0.0,1.0,0.01,1.0e-06,X,3,M,F)
WRITE(*,10) M
10 FORMAT(1X,'M=',I2)
DO 20 I=1,M
20 WRITE(*,30) I,X(I)
30 FORMAT(1X,'X(',I2,')=',E25.6)
END
FUNCTION F(X)
F=-2500/(1+x)-2500/(1+x)**2+830/(1+x)**3+1800/(1+x)**4+1600/(1+x)
$**5+1900*((1+x)**4-1)/(x*(1+x)**4*(1+x)**5)+3840/(1+x)**10
RETURN
END

编译可以通过，但是运行时出现“浮点错误”。