      program AM
      implicit real*8 (a-h,k,o-z)
      parameter (N=10)
      dimension x(0:N),y(0:N)
      
      x(0)= 0d0
      x(N)= 1d0
      h=(x(N)-x(0))/N
      print*,'N   =',N
      print*,'x(0)=',x(0)
      print*,'x(N)=',x(N)
      print*,'h   =',h
      print*
      
      do i=0,N
         x(i)=i*h+x(0)
      enddo
      
      y(0)=1d0
      print '(3x,''i'',14x,''x'',14x,''y'',11x,''y_an'')'
*     For the first steps we use a single-step method, e.g. the 4th-order Runge-Kutta method
      do i=0,1
         k1=h*f(x(i)    ,y(i))
         k2=h*f(x(i)+h/2,y(i)+k1/2)
         k3=h*f(x(i)+h/2,y(i)+k2/2)
         k4=h*f(x(i)+h  ,y(i)+k3  )
         y(i+1)=y(i)+(k1+2*k2+2*k3+k4)/6
      enddo
*     Then, we use the multistep method
      do i=2,N-1
         y(i+1)=y(i)+h/12*(23*f(x(i),y(i))-16*f(x(i-1),y(i-1))
     &                                     +5*f(x(i-2),y(i-2)))
      enddo
      
      do i=0,N
         print '(i4,3f15.10)',i,x(i),y(i),2d0*dexp(x(i))-x(i)-1d0
      enddo
      
      pause
      end
      
      function f(x,y)
      implicit real*8 (a-h,o-z)
      f=x+y
      return
      end
