***LEGENDRE POLYNOMIAL***********************************************
*      PROGRAM FOR GENERATING LEGENDRE POLYNOMIAL OF ANY DEGREE     *
*********************************************************************
      program legendrepoly
      implicit doubleprecision (a-h,o-z)
*
      print*, 'Input the degree of polynomial:'
      read(*,*)n  
*     n = 4
      print*, 'Input the starting value of independent variable, x:'
      read(*,*)xl
*     xl = -1
      print*, 'Input the ending value of independent variable, x:'
      read(*,*)xf
*     xf = 1
      print*, 'Input the print interval delta x:'
      read(*,*)dx
*     dx = 1d-1
      x = xl
      write(*,*)
      write(*,101)'x','P(',n,')'
10    if ( x <= (xf+(dx/2)) ) then
         call bonnet(n,x,p)
         write(*,102)x,p
         x = x + dx
         goto 10
      endif
*
101   format(a8,a13,i1,a1)
102   format(1x,e12.5,1x,e14.7)
      stop
      end
*
***********************************************************************
*
      subroutine bonnet(n,xx,p)
      implicit doubleprecision (a-h,o-z)
      dimension pleg(0:10)
c
c     Legendre polynomial
      pleg(0) = 1d0
      pleg(1) = xx
      do i = 2,n
         pleg(i) = ( (2*i-1)*xx*pleg(i-1) - (i-1)*pleg(i-2) ) /i
      enddo
      p = pleg(n)
c
      return
      end