program newtfind
      
c   program to find root of Rayleigh wave cubic

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

      print *,' input initial values of damp,x'
5     read (5,*,end=90) damp,x0
      x=x0
      do 10 i=1,100
        j=i
        call quartic(x,f,df)
        x1=x-damp*f/df
        print *,x1, ' new estimate of root'
        print *,'enter new damping'
        read (5,*) damp
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)
      implicit none
      real x,f,df
      f=3*x**4-16*x**3+18*x*x 
      df = 12*x*(x-1)*(x-3)
      return
      end