SUBROUTINE ORTHFE(X,P,Y,DY)
*     Return function value Y and error DY of the orthogonal polynomial
*     at the abscissa value X, using the array P (see ORTHFT)                 
      REAL X,Y,DY,P(5,*)
      DOUBLE PRECISION Q(2),F,DF
      NP=MOD(P(1,1),100.0)+0.5
      Q(1)=P(3,1)       ! constant term
      Q(2)=0.0
      F =P(3,1)*P(4,1)  
      DF=P(3,1)*P(3,1)
      IAM=1
      DO K=2,NP
       IAM=3-IAM  ! recurrence relation for higher terms
       Q(IAM)=((X-P(1,K))*Q(3-IAM)-P(2,K)*Q(IAM))*P(3,K)
       F =F +P(4,K)*Q(IAM) 
       DF=DF+Q(IAM)*Q(IAM)
      END DO 
      Y =F               ! function value
      DY=SQRT(DF*P(2,1)) ! standard deviation 
      END