Program 1 : fd1, fdm, spd

      program prog1 
      implicit none

      real cubic,slope,rms
      external cubic,slope,rms  
      real f(10),x(10),h(10),fslp(10),sslp(10),sslp1,sslpn
      real a(10),c(10),d(10)
      integer i

# define vector of arbitrary coordinates, x,
# compute values of f and slope of f at these coordinates
# compute spacing of coordinates
      print *,'The cubic f=x**3+15*x**2+60*x+6:'
      do i=1,10
         x(i)=i**2
         f(i)=cubic(x(i))
         print *,'f(x(',i,' ))=',f(i),'    x(',i,' )=',x(i)
         fslp(i)=slope(x(i))
         if (i.ne.1) h(i-1)=x(i)-x(i-1)
      enddo
      
      call fd1(sslp1,f,h)   
      print *,''
      print *,'The slope at the first point'
      print *,'  computed by fd1:',sslp1
      print *,'  computed analytically:',fslp(1)

      call fdm(sslpn,f,h,10)
      print *,''
      print *,'The slope at the last point'
      print *,'  computed by fdm:',sslpn
      print *,'  computed analytically:',fslp(10)
      print *,''

      call spd(sslp,f,10,h,0,sslp1,0,sslpn,a,c,d)
      print *,'For the ten slopes calculated,'
      print *,'spd has an rms error of:',rms(sslp,fslp,10)
      
      stop
      end
   
      real function cubic(x)
      real x
      cubic=x**3+15*x**2+60*x+6
      return
      end
         
      real function slope(x)
      real x
      slope=3*x**2+30*x+60
      return
      end

      real function rms(a,b,n)
      real a(n),b(n),sum
      integer i,n
      sum=0
      do i=1,n
         sum=sum+(a(i)-b(i))**2
      enddo
      rms=sqrt(sum/n)
      return
      end