program prog05
      implicit real*8 (a-h,o-z)
c
c     This program solves the differential equation
c      d^2 x/d t^2 = -g sin(x)
c     in a improved Euler method
c 
c Set default parameters
      ti = 0.0d0         ! Initial Value of t
      tf = 10.0d0        ! Final   Value of t
      xi = 2.0d0         ! Inital  Value of x
      vi = 0.0d0         ! Inital  Value of v
c Parameter inputs
      write(*,*) '# n=?'
      read(*,*) n
c
      dt = (tf-ti)/n     ! Time Mesh
      t=ti               ! Inital Condition
      x=xi               ! Inital Condition
      v=vi               ! Inital Condition
c
         write( *,800) t,x,v,energy(x,v)
      do i = 1,n
         v0 = v                ! dx/dt at (x,v)
         f0 = f(x)             ! dv/dt at (x,v)
c
         x1 = x + dt*v0
         v1 = v + dt*f0        ! dx/dt at (x1,v1) 
         f1 = f(x1)            ! dv/dt at (x1,v1) 
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
 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