c***********************************************************************
c                                                                      *
c we use the subroutine dverk (with default options) to solve the      *
c simple test problem                                                  *
c      y'(1) = -y(2),  y(1) = 1 at x = 0,                              *
c      y'(2) =  y(1),  y(2) = 0 at x = 0.                              *
c this test program prints the computed solution  and  the  associated *
c error at the output points x = 1, 2, 3, ... 10.                      *
c                                                                      *
c this program contains some fortran 77 features.                      *
c k.r. jackson   11/14/85                                              *
c                                                                      *
c***********************************************************************
c
        integer n, nw, ind, k
        double precision x, xend, y(2), tol, c(24), w(2,9)
        double precision error(2)
        external fcn
c
c initialize the variables.
        n   = 2
        nw  = 2
        ind = 1
c
        y(1) = 1.d0
        y(2) = 0.d0
        x    = 0.d0
c
        tol = 1.d-6
c
c output the heading.
        print 1 
 1      format(/// ' output from simple test program for dverk.'
     +     // 6x, 'x', 6x, 'y(1)', 11x, 'error(1)', 
     +                 7x, 'y(2)', 11x, 'error(2)' 
     +      / 5x, '---', 4(3x, '------------') )
c
c calculate and output the solution and the associated error
c at the output points x = 1, 2, 3, ... 10.
c
        do 10 k = 1, 10
           xend = dfloat(k)
           call dverk ( n, fcn, x, y, xend, tol, ind, c, nw, w )
           if ( ind .eq. 3 ) then
c             integration successful - compute errors and print results
              error(1) = dcos(x) - y(1)
              error(2) = dsin(x) - y(2)
              print 2, x, y(1), error(1), y(2), error(2)
 2            format( 4x, f4.1, 4d15.6 )
           else
c             integration unsuccessful - print message and stop
              print 3, x, ind
 3            format(/// ' integration stopped at x =', f7.3,
     +                   ' with ind =', i5 )
              stop
           end if
 10     continue
c
        stop
        end
c
c
c
        subroutine fcn ( n, x, y, yprime )
c
        integer n
        double precision x, y(n), yprime(n)
c
        yprime(1) = -y(2)
        yprime(2) =  y(1)
c
        return
        end