program ex053
      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='ex053.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
c...
c Select one of the following methods.
         call impeul(x,v,t,dt,gam,omg)
c...
         write( *,800) t,x,v
         write(16,800) t,x,v
      end do
c
      stop
 800  format(3(1x,f15.7))
      end
c *********************************************************************
      subroutine impeul(x,v,t,dt,gam,omg)
c *********************************************************************
      implicit real*8 (a-h,o-z)
c
      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
      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