c
c    Comparisons of various IVP solvers
c       Euler, Backward Euler, Leapfrog, Trapezoidal
c
c    IVP:  y' = r*(1-y)  where  y(0)=a
c
	implicit real*8(a-h,o-z)
      dimension t(0:10000), y(0:10000), e(0:10000), be(0:10000)
	dimension frog(0:10000),tt(0:10000),trap(0:10000)
	parameter(r=6.0)
c      f(t,y) = 1 + r*y
c      exact(t)= - (1.0-exp(r*t))/r
	f(t,y)=r*(1-y)
	exact(t)=1-exp(-r*t)
	do 80 in=1,8
	n=2**(in+2)
      a=0
      tmax=1.0
      t(0)=0
	tt(0)=0
      y(0)=a
      e(0)=a
      be(0)=a
	frog(0)=a
	trap(0)=a
      dt=tmax/n
*
* Euler
*
      do 10 j=1,n
        t(j)=dt*j
   10   e(j)=e(j-1)+dt*f(t(j-1),e(j-1))
*
* Backward Euler and Trapezoidal
*
      do 20 j=1,n
        q0=be(j-1)
        do 16 jj=1,100
          q=be(j-1)+dt*f(t(j),q0)
          if(abs((q-q0)/q).lt.0.00001) go to 18
   16     q0=q
   18	  trap(j)=0.5*(q+e(j))
   20   be(j)=q
*
* Leapfrog
*
	frog(1)=0.5*(e(1)+be(1))
	do 40 j=2,n
   40	  frog(j)=frog(j-2)+2*dt*f(t(j-1),frog(j-1))
*
*  Exact
*
	h=tmax/1000
	do 60 j=1,1000
	  tt(j)=h*j
   60	  y(j)=exact(tt(j))
*
* Output
*
c      write(6,100) (t(j),e(j),be(j),frog(j),trap(j),j=0,n)
  100 format(f10.3,',',e12.4,',',e12.4,',',e12.4,',',e12.4)
	write(6,101) (tt(j),y(j),j=0,1000)
  101	format(f10.3,',',e12.4)
*
* Error
*
	error1=abs(y(1000)-e(n))
	error2=abs(y(1000)-be(n))
	error3=abs(y(1000)-frog(n))
	error4=abs(y(1000)-trap(n))
c	write(6,103) n,error1,error2,error3,error4
  103 format(i5,',',e12.4,',',e12.4,',',e12.4,',',e12.4)
   80 continue
      end