program sine
cc
cc  Program written by J. M. Powers
cc  University of Notre Dame
cc  1999
cc
cc  This program solves the second order ODE:
cc
cc           d2f/dx2  +  f  = 0
cc
cc  subject to the initial conditions
cc
cc           f(0) = 0,  df/dx(0) = 1.
cc
cc  The solution is known to be 
cc
cc              f = sin(x).
cc
cc  The program compares its estimation to the known solution.
cc  The program calls upon the IMSL subroutine DIVPRK
cc  which uses a Runge-Kutta-Verner fifth and sixth order
cc  explicit method to numerically integrate the resulting
cc  coupled set of first order ODE's.
cc
cc  IMSL documenation is available with the SUN command
cc
cc     iptdoc
cc
cc  On the Sun system at Notre Dame, the proper compile
cc  statement is
cc
cc    f77 sin.f -o sin -limsl -lnsl -lsocket
cc
cc  This generates an executable file named "sin"
cc
cc
cc  It's a good idea to declare "implicit none", which
cc  forces you to declare all variables.
cc
      implicit none
cc
      double precision fcn,x,xend,tol,param(50),f(2)
      double precision dx,xmax
      integer ido,neqn,i
      external fcn
cc
cc  Note to be in double precision, all variables must be
cc  declared double precision; AND all constants must
cc  have a "d0" extension.

      dx = 0.01d0
      xmax = 5.d1
      tol = 1.d-10
      x = 0.d0
      param(4) = 10**8
      ido = 1
      neqn = 2
cc
cc  Initial conditions
cc
      f(1) = 0.d0
      f(2) = 1.d0
cc
cc  Print headers

cc
      print110,'x','f(x)','g(x)','sin(x)'
110   format(a7,3a15)
      print*,' '
cc
cc Print initial values
cc
      print100,x,f(1),f(2),dsin(x)
100   format(4(e12.6,3x))

cc
cc  Main loop in time
cc

      do 10 i=1,100000
	 xend = x + dx
cc
cc  Below is the IMSL routine.
cc
	 call divprk(ido,neqn,fcn,x,xend,tol,param,f)
	 x = xend
	 print100,x,f(1),f(2),1.d0-f(2)**2-f(1)**2
	 if(x.ge.xmax)go to 11
10       continue

11    continue

      stop
      end

cc
cc  Subroutine fcn below is needed by
cc  the IMSL routine divprk to calculate derivatives.
cc
      subroutine fcn(neqn,x,f,fprime)

      implicit none
      double precision x,f(2),fprime(2)
      integer neqn

      fprime(1) =  f(2)
      fprime(2) = -f(1)

      return
      end