      program ex052
      implicit real*8 (a-h,o-z)
c 
c Set default parameters
      n = 1000           ! mesh points
      gam = 0.1d0        ! gamma, Friction (Masatsu)
      omg  = 2.0d0       ! omega, Restoring Force
      ti = 0.0d0         ! Initial Value of t
      tf = 100.0d0       ! Final   Value of t
      xi = 5.0d0         ! Initial  Value of x
      vi = 0.0d0         ! Initial  Value of v
C
c Parameter inputs
      write(*,810) '#  gam, omg =?'
      read(*,*)        gam, omg
 810  format(a)
      open(16,file='ex052.dat') ! open(unit,file='filename')
c
      dt = (tf-ti)/n     ! Time Mesh
      t=ti               ! Initial Condition
      x=xi               ! Initial Condition
      v=vi               ! Initial Condition
c
      write( *,800) t,x,v
      write(16,800) t,x,v
      do i = 1,n
         v0 = v
         f0 = f(v,x,t,gam,omg)
c
         t1 = t + dt
         x1 = x + dt*v0
         v1 = v + dt*f0
         f1 = f(v1,x1,t1,gam,omg)
c
         dx = ( v0 + v1 )*dt/2	! Improved Euler
         dv = ( f0 + f1 )*dt/2	! Improved Euler
c
         t = t + dt
         x = x + dx
         v = v + dv
c
         write( *,800) t,x,v
         write(16,800) t,x,v
      end do
c
      stop
 800  format(3(1x,f15.7))
      end
c **********************************************************************
      function f(v,x,t,gam,omg)
c **********************************************************************
      implicit real*8 (a-h,o-z)
c Force
      f = -gam*v-omg*omg*x+cos(t)
      end
