c*****************************************************************

	SUBROUTINE maehly(a,n,nr,ki)

c recherche des zéros réels par l'algorithme de Maehly

c entrées:
c	a(i): coefficient de x**i du polynôme de degré <= 30
c	n: degré
c	nr<=n : nombre de racines à obtenir

c sorties
c	ki:racines

c Auteur P. Morel laboratoire Cassiopée, OCA, CESAM2k

c---------------------------------------------------------------------

	USE mod_kind

	IMPLICIT NONE

	REAL(kind=dp), INTENT(in), DIMENSION(0:n) :: a
	REAL(kind=dp), INTENT(out), DIMENSION(n) :: ki
	INTEGER, INTENT(in) :: n, nr

	REAL(kind=dp), DIMENSION(0:n) :: p

	REAL(kind=dp), PARAMETER :: epsi=1.d-12
	REAL(kind=dp), SAVE :: x0
	REAL(kind=dp) :: eps, s, x

	INTEGER, PARAMETER  ::  der=1, nmax=100
	INTEGER :: compt, i, j, k

	LOGICAL, SAVE :: init=.TRUE.

c--------------------------------------------------------------------------

2000	FORMAT(8es10.3)

c initialisation
	IF(init)THEN
	 init=.FALSE.
	 x0=EXP(1.d0)
	ENDIF

	DO j=1,nr
	 x=x0
	 eps=1.d0
	 compt=0
	 B1: DO
	  IF(ABS(eps) <= epsi .OR. compt > nmax)EXIT B1
	  compt=compt+1 ; s=0.d0
	  DO i=1,j-1
	   s=s+1.d0/(x-ki(i))
	  ENDDO
	  CALL polyder(a,n,der,x,p)
	  eps=p(0)/(p(1)-p(0)*s) ;  x=x-eps
	 ENDDO B1

c sortie
	 IF(compt >= nmax)THEN
	  PRINT*,nmax,' iterations et pas de convergence dans maehly'
	  PRINT*,'Nombre de racines deja obtenues',j-1
	  DO k=1,j-1
	   PRINT*,k,ki(k)
	  ENDDO
	  STOP
	 ELSE
	  ki(j)=x
	 ENDIF
	ENDDO

	RETURN

	END SUBROUTINE maehly