subroutine approx(x,y,yd,kdim)
c
c     The routine performs third-order interpolation between points
c     yold and ynew parameterized by an independent variable x.
c     dold, dnew are the first derivatives of y with respect
c     to x at the points yold and ynew.
c
      dimension y(18),yd(18)
      common/appr/ xold,xnew,yold(18),dold(18),ynew(18),dnew(18)
c
      a=(x-xnew)/(xnew-xold)
      aux=a+1.
      a1=(2.*a+3.)*a*a
      a2=1.-a1
      b1=aux*a*(x-xnew)
      b2=aux*a*(x-xold)
      ad1=6.*a*aux/(xnew-xold)
      ad2=-ad1
      bd1=a*(3.*a+2.)
      bd2=aux*(3.*a+1.)
      do 1 i=1,kdim
      y(i)=a1*yold(i)+a2*ynew(i)+b1*dold(i)+b2*dnew(i)
      yd(i)=ad1*yold(i)+ad2*ynew(i)+bd1*dold(i)+bd2*dnew(i)
    1 continue
      return
      end