program test
      parameter( n=10 )
      parameter( nmax=100 )
      real p(0:n,nmax)
      real x(nmax), y(nmax), theta(nmax)
      real tmin, tmax
      real dt, pi
      pi=3.14159
      tmin=-pi/2.0
      tmax=pi/2.0
      dt=(tmax-tmin)/(nmax-1)
      write(*,*) "input the number, show number"
      read(*,*) k,l
      do 10 i=1,nmax
         theta(i)=tmin+dt*(i-1)
         x(i)=cos(theta(i))
  10  continue
      do 30 i=1,n
         do 40 j=1,nmax
            p(i,j)=0.0
  40  continue
  30  continue
      CALL LEGENDRE(k,nmax,x,p(0:k,1:nmax))
      do 20 i=1,nmax
         y(i)=p(l,i)
  20  continue

*-- グラフ ----
      WRITE(*,*) ' WORKSTATION ID (I)  ? ;'
      CALL SGPWSN
      READ (*,*) IWS

      CALL GROPN( IWS )
      CALL GRFRM

      CALL USSPNT( nmax, theta, y )
      CALL GRSWND( tmin, tmax, -1.0, 1.0 )
      CALL GRSVPT(  0.2,  0.8,  0.2,  0.8 )
      CALL USPFIT
      CALL GRSTRF

      CALL USDAXS
      CALL UULIN( nmax, theta, y )

      CALL GRCLS

      stop
      end

*-- ルジャンドル漸化式のサブルーチン ---
      subroutine LEGENDRE(n, nmax, x, p)
      real p(0:n,nmax), x(nmax)
*-- 初項の設定 ---
      do 10 i=1,nmax
         p(0,i)=1.0
  10  continue
      if(n.eq.0)then
         go to 100
      end if
      do 20 i=1,nmax
         p(1,i)=x(i)
  20  continue
      if(n.eq.1)then
         go to 100
      end if
*-- 漸化式の計算 ---
      do 30 j=1,n-1
         do 40 i=1,nmax
            p(j+1,i)=(1.0/real(j+1))*(p(j,i)*(2.0*j+1)*(p(1,i))
     &               -real(j)*p(j-1,i))
  40  continue
  30  continue
  100 return
      end