       program pen_pla
       implicit logical (a-z)
       double precision theta, omega,pi
       double precision y(2),dy(2), omega0, theta0, delta, time, omega00
       integer i, j, k
       pi=4.d0*datan(1.d0)
       omega00=5.d-1
       delta=1.d-3
       do k=-2,2
       omega0=dble(k)*omega00
       do j=1,21,5
         theta0=-pi+dble(j-1)*pi/10.d0
         time=0.d0
         y(1)=theta0
         y(2)=omega0
         do i = 1,15000
           call eco(y,dy)
           write(*,*) time, y(1), y(2)
           call rk(y,dy,time,delta,y)
           time=time+delta
         end do
         write(*,*) ' '
       end do
         write(*,*) ' '
       end do
       end program
      subroutine eco(y,dy)
      implicit logical (a-z)
      double precision omega,theta, y(2), dy(2)
      theta=y(1)
      omega=y(2)
      dy(1)=omega
      dy(2)=-sin(theta)
      return
      end  
      subroutine rk(y,dydx,x,delta,yout)
      implicit logical (a-z)
      integer i
      double precision y(2),dydx(2),yout(2),yt(2),dyt(2),dym(2)
      double precision hh,delta,h6,x,xh
      hh=delta*0.5d0
      h6=delta/6.d0
      xh=x+hh
      do i=1,2
        yt(i)=y(i)+hh*dydx(i)
      end do
      call eco(yt,dyt)
      do i=1,2
        yt(i)=y(i)+hh*dyt(i)
      end do
      call eco(yt,dym)
      do i=1,2
        yt(i)=y(i)+delta*dym(i)
        dym(i)=dyt(i)+dym(i)
      end do
      call eco(yt,dyt)
      do i=1,2
        yout(i)=y(i)+h6*(dydx(i)+dyt(i)+2.d0*dym(i))
      end do
      return
      end
