chapter ii. runge example
c  from  * a practical guide to splines *  by c. de boor    
      integer i,istep,j,k,n,nmk,nm1
      real aloger,algerp,d(20),decay,dx,errmax,g,h,pnatx,step,tau(20),x
      data step, istep /20., 20/
      g(x) = 1./(1.+(5.*x)**2)
      print 600
  600 format(28h  n   max.error   decay exp.//)
      decay = 0.
      do 40 n=2,20,2
c        choose interpolation points  tau(1), ..., tau(n) , equally
c        spaced in (-1,1), and set  d(i) = g(tau(i)), i=1,...,n.
         nm1 = n-1
         h = 2./float(nm1)
         do 10 i=1,n
            tau(i) = float(i-1)*h - 1.
   10       d(i) = g(tau(i))
c        calculate the divided differences for the newton form.
c
         do 20 k=1,nm1
            nmk = n-k
            do 20 i=1,nmk
   20          d(i) = (d(i+1)-d(i))/(tau(i+k)-tau(i))
c
c        estimate max.interpolation error on (-1,1).
         errmax = 0.
         do 30 i=2,n
            dx = (tau(i)-tau(i-1))/step
            do 30 j=1,istep
               x = tau(i-1) + float(j)*dx
c              evaluate interp.pol. by nested multiplication
c
               pnatx = d(1)
               do 29 k=2,n
   29             pnatx = d(k) + (x-tau(k))*pnatx
c
   30          errmax = amax1(errmax,abs(g(x)-pnatx))
         aloger = alog(errmax)
         if (n .gt. 2)  decay =
     *       (aloger - algerp)/alog(float(n)/float(n-2))
         algerp = aloger
   40    print 640,n,errmax,decay
  640 format(i3,e12.4,f11.2)
                                        stop
      end