program newtfind1
      
c   program to find root of Rayleigh wave cubic

      implicit none
      integer i,j
      real x0,x,x1,f,df,ddf
      character*1 yn

      print *,' input initial value of x'
5     read (5,*,end=90) x0
      x=x0
      do 10 i=1,100
        j=i
        call quartic(x,f,df,ddf)
        print *,x,f,df,ddf
        x1=x-df/ddf
        print *,x1, ' new estimate of root'
c       print *, ' stop yet?'
c       read (5,*) yn
c       if (yn .eq. 'y') stop
        if (abs(x-x1) .lt. 1.e-4) then
          print *,j,x1 
          go to 5
        end if
        x=x1
10    continue
        go to 5
90      stop
        end

      subroutine quartic(x,f,df,ddf)
      implicit none
      real x,f,df,ddf
      f=3*x**4-16*x**3+18*x*x 
      df = 12*x*(x-1)*(x-3)
      ddf=36*x*x-96*x+36
      return
      end