parameter (N=64, NN=2*N+1)
       implicit real*8(a-h,o-z)
       real x(0:NN+1),ubar(0:NN+1,0:3)
       real uapp1(0:NN+1), xL(0:NN+1,0:3)
       complex abar(0:NN+2), al(0:NN+2), si(0:NN+1)
       complex w, A, y, ss, ss1, ss2, sum, sum1

       pi=4.0*atan(1.0)
       w=cmplx(0, 1)
       alpha=15.0*log(10.0)

       dx=2.0*pi/NN
       dx1=0.5*dx
       dt=0.01
       ttime=0.0
       Tmax=0.5

       par=1.9

c compute the initial fourier coefficient
       do j=0, NN
          x(j)=j*dx
          if ( (x(j)-dx1) .lt. par .and. (x(j)+dx1) .gt. par) jp=j
       enddo

       do j=0, NN
          x(j)=j*dx
c          ubar(j, 0)=0.3+0.7*(-cos(x(j)+0.5*dx)+cos(x(j)-0.5*dx))/dx

          s1=cos(0.5*(x(j)+dx1))
          s2=cos(0.5*(x(j)-dx1))

          if (j .lt. jp) then
             ubar(j, 0)=-2.0*(s1-s2)/dx
           elseif (j .gt. jp) then
             ubar(j, 0)=2.0*(s1-s2)/dx
           else
             p1=-2.0*(cos(0.5*par)-s2)/dx
             p2=2.0*(s1-cos(0.5*par))/dx
             ubar(j, 0)=p1+p2
          endif

       enddo

 99    do NR=0, 2

       do k=0, NN-1
           k1=k-0.5*(NN-1)
           ss=0.0
           do j=0, NN-1
              s1=cos(k1*x(j))
              s2=-sin(k1*x(j))
              ss=ss+ubar(j,NR)*cmplx(s1,s2)
           enddo
           abar(k)=ss/NN
       enddo

       do k=0, NN-1
           k1=k-0.5*(NN-1)
           if (k1 .eq. 0) then
           si(k)=1.0
           else
           si(k)=sin(0.5*k1*dx)/(0.5*k1*dx)
           endif
           al(k)=abar(k)/si(k)
       enddo

      do j=0, NN-1
          xj=(j+0.5)*dx
          sum=0.0

          do k=1, N
             tau=exp(-alpha*(1.0*k/N)**20)

             ps=k*xj
             s1=cos(ps)
             s2=sin(ps)
             ss2=cmplx(s1, s2)
             ss=al(k+N)*ss2
             sum=sum+tau*ss

             k1=-k
             ps=k1*xj
             s1=cos(ps)
             s2=sin(ps)
             ss2=cmplx(s1, s2)
             ss=al(k1+N)*ss2

             sum=sum+tau*ss
          enddo

          uapp1(j)=real(al(N)+sum) 

       enddo

       do j=0, NN-1
          ss1=uapp1(j)
          if (j .eq. 0) then
             ss2=uapp1(NN-1)
             else
             ss2=uapp1(j-1)
          endif
          xL(j, NR)=-0.5*(ss1*ss1-ss2*ss2)/dx

          if (NR .eq. 0) then
             ubar(j,1)=ubar(j,0)+dt*xL(j,0)
          endif

          if (NR .eq. 1) then
             sa=0.75*ubar(j,0)+0.25*ubar(j,1)+0.25*dt*xL(j,1)
             ubar(j,2)=sa
          endif

          if (NR .eq. 2) then
             sa=ubar(j,0)/3.0+ubar(j,2)*2/3.0+dt*xL(j,2)*2/3.0
             ubar(j,3)=sa
          endif

        enddo

       enddo

       ttime=ttime+dt
       print *, ttime, '    time'
       if (ttime .lt. Tmax) then
          do j=0, NN
            ubar(j,0)=ubar(j, 3)
          enddo
          goto 99
       endif

       do k=0, NN-1
           k1=k-0.5*(NN-1)
           ss=0.0
           do j=0, NN-1
              p1=k1*x(j)
              s1=cos(p1)
              s2=-sin(p1)
              ss=ss+ubar(j,NR)*cmplx(s1,s2)
           enddo
           abar(k)=ss/NN
           al(k)=abar(k)/si(k)
       enddo

      do j=0, NN-1
          xj=x(j)
          sum=0.0

          do k=1, N
             tau=exp(-alpha*(1.0*k/N)**4)

             ps=k*xj
             s1=cos(ps)
             s2=sin(ps)
             ss2=cmplx(s1, s2)
             ss=al(k+N)*ss2
             sum=sum+tau*ss

             k1=-k
             ps=k1*xj
             s1=cos(ps)
             s2=sin(ps)
             ss2=cmplx(s1, s2)
             ss=al(k1+N)*ss2

             sum=sum+tau*ss
          enddo

          uapp1(j)=real(al(N)+sum)

       enddo

       write (17, *) 'x=['
       do j=0, NN-1
          write (17, *) x(j)
       enddo
       write (17, *) '];'

        write (18, *) 'y=['
       do j=0, NN-1
          write (18, *) uapp1(j)
       enddo
       write (18, *) '];'

       print *, A, y, '   A, y'
       
       stop
       end