c===========================================================
c     Self-interaction potentials routine
c===========================================================
      subroutine pot(neq,y,u,dudphi2)
         implicit    none

         include    'fcn.inc'
         integer     neq
         integer     i
         real*8      y(neq)

         real*8      u,     dudphi2




      if (pty.eq.0) then  ! massless KG field

           u = 0.0d0
      dudphi2 = 0.0d0
         
      else if (pty.eq.1) then  !polynomial self-interacting BS 

      if(pdim.le.100)then
            u = 0.0d0
      dudphi2 = 0.0d0 
            do i = 1, pdim
                    u = u + pvec(i)*y(3)**i   ! y(3) = phi0(r) = |phi(r)|
              dudphi2 = dudphi2 + 0.5d0*i*pvec(i)*y(3)**(i-2)
            enddo
c            write(*,*)y(3), dudphi2 
c            write(*,*)y(3), u
      else
        write(0,*)' pot: The maximal value allowed for pn, dimension'
        write(0,*)'      of the parameter vector p, is 100'
        write(0,*)' '
        write(0,*)' To modify the previous value, change the common'
        write(0,*)' statement at fcn.inc and the if statement at pot.f'
        write(0,*)' For example:'
        write(0,*)' '
        write(0,*)'              p(100)   ------->   p(1000)'
        write(0,*)'  if(pdim.le.100)then -------> if(pdim.le.1000)then'
      endif

              

      else

      write(*,*)' Wrong selection of the potential!'
      write(*,*)' pt should be an integer in the range [0,1]'
      write(*,*)
      stop

      endif
            
         return
      end