program Cheby
      parameter( n=50 )
      parameter( nmax=200 )
      real che(0:n,nmax)
      real x(nmax), y(nmax), theta(nmax)
      real tmin, tmax
      real dx, pi, dt
      pi=3.14159
      tmin=-pi
      tmax=pi
      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
            che(i,j)=0.0
  40  continue
  30  continue
      CALL CHEBYSHEV(k,nmax,che(0:k,1:nmax))
      do 20 i=1,nmax
         y(i)=che(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 CHEBYSHEV(n, nmax, che)
      real che(0:n,nmax)
      real xmin, xmax
      real dx
      xmin=-1.0
      xmax=1.0
      dx=(xmax-xmin)/(nmax-1)
*-- 初項の設定 ---
      do 10 i=1,nmax
         che(0,i)=1.0
  10  continue
      if(n.eq.0)then
         go to 100
      end if
      do 20 i=1,nmax
         che(1,i)=xmin+dx*(real(i)-1.0)
  20  continue
      if(n.eq.1)then
         go to 100
      end if
*-- 漸化式の計算 ---
      do 30 j=1,n-1
         do 40 i=1,nmax
            che(j+1,i)=2.0*che(1,i)*che(j,i)-che(j-1,i)
  40  continue
  30  continue
  100 return
      end