      program ex051
      implicit real*8 (a-h,o-z)
c 
c Set default parameters
      n = 1000           ! mesh points
      ti = 0.0d0         ! Initial  Value of t
      tf = 10.0d0        ! Final    Value of t
      xi = 0.5d0         ! Initial  Value of x
      vi = 0.0d0         ! Initial  Value of v

c Parameter inputs
      write(*,810) '# n, xi, vi =?'
      read(*,*) n, xi, vi
 810  format(a)
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,energy(x,v)
      do i = 1,n
         v0 = v
         f0 = f(x)
c
         x1 = x + dt*v0
         v1 = v + dt*f0
         f1 = f(x1)
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
         write( *,800) t,x,v,energy(x,v)
      end do
c
      stop
 800  format(4(1x,f15.7))
      end
c **********************************************************************
      function f(x)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      g = 9.8d0          ! gravitation constant
      f=-g*sin(x)        ! Force
      end
c **********************************************************************
      function energy(x,v)
c **********************************************************************
      implicit real*8 (a-h,o-z)
      g = 9.8d0          ! gravitation constant
      energy = v*v/2 + g*(1.0d0-cos(x))  ! Energy
      end
