c @(#)Copyright (c), 1987, 1993 StatSci, Inc.  All rights reserved.
C       Heun's method for solving dy/dt=f(t,y), using step size h :
C       k1 = h f(t,y)
C       k2 = h f(t+h,y+k1)
C       ynext = y + (k1+k2)/2

        subroutine heun(f, neqn, y, tstart, tend, step, work)
        integer neqn
        real*4 f, y(neqn), tstart, tend, step, work(neqn,3)
C       work(1,1) is k1, work(1,2) is k2, work(1,3) is y+k1
        integer i, nstep, istep
        real*4 t
        external f

        nstep = max((tend - tstart) / step, 1.0)
        step = (tend - tstart) / nstep
        do 30 istep = 1, nstep
                t = tstart + (istep-1) * step
                call f(t, y, work(1,1))
                do 10 i = 1, neqn
                        work(i,1) = step * work(i,1)
                        work(i,3) = y(i) + work(i,1)
10              continue
                call f(t+step, work(1,3), work(1,2))
                do 20 i = 1, neqn
                        work(i,2) = step * work(i,2)
                        y(i) = y(i) + 0.5 * (work(i,1) + work(i,2))
20              continue
30      continue
        return
        end