c   demonstration program for ode/de/step package.
c
c   reference:shampine and gordon's computer solution of
c     ordinary differential equations:the initial value problem.
c
c   the package is used to solve the defining equations for the
c   jacobian elliptic functions:
c         y1'=y2*y3,      y1(0)=0
c         y2'=-y1*y3,     y2(0)=1
c         y3'=-ksq*y1*y2, y3(0)=1
c   which is solved for ti=i*5 for i=0, 1, ..., 12.
c   the analytic solution is
c         y1=sn(t/ksq)
c         y2=cn(t/ksq)
c         y3=dn(t/ksq)
c
c   in this case we use ksq=k*k=0.51.
c   output is written on logical unit lout which has been set to 6.
c
      real y, t, tout, relerr, abserr, work
      dimension y(3), work(163), iwork(5)
      external f
      data lout/6/
c
      neqn=3
      relerr=1.0e-9
      abserr=1.0e-16
      iflag=1
      write(lout, 1000)neqn, relerr, abserr, iflag
c
      t=0.0e0
      y(1)=0.0e0
      y(2)=1.0e0
      y(3)=1.0e0
      write(lout, 1002)t, y, iflag
      do 11 i=1,12
      tout=5.0e0*i
    3 call ode(f, neqn, y, t, tout, relerr, abserr, iflag, work,
     * iwork)
      write(lout, 1002)t, y, iflag
      goto (900,11,5,7,9,900),iflag
    5 write(lout, 1004)
      goto 3
    7 write(lout, 1006)
      goto 900
    9 write(lout, 1008)
      goto 900
   11 continue
  900 stop
 1000 format('- neqn=',i3,'  relerr=',1pe12.2,'    abserr=',f12.2,
     * '   iflag=',i3,//)
 1002 format(' ',4(1pe16.8,2x),i3)
 1004 format('0tolerances too much and have been changed')
 1006 format('0too many steps...eqn. is hard')
 1008 format('0eqn. appears to be stiff')
      end
      subroutine f(t, y, yp)
      dimension y(3), yp(3)
      real t, y, yp
      yp(1)=y(2)*y(3)
      yp(2)=-y(1)*y(3)
      yp(3)=-0.51e0*y(1)*y(2)
      return
      end