function sum=gaussint7(f,a,b)
%
%              function sum=gaussint7(f,a,b)
%
%           the abscissae and weights are given for the interval (-1,1).
%           because of symmetry only the positive abscissae and their
%           corresponding weights are given.
%
%           xgk    - abscissae of the 15-point kronrod rule
%                    xgk(2), xgk(4), ...  abscissae of the 7-point
%                    gauss rule
%                    xgk(1), xgk(3), ...  abscissae which are optimally
%                    added to the 7-point gauss rule
%
%           wgk    - weights of the 15-point kronrod rule
%
%           wg     - weights of the 7-point gauss rule
%
%
% gauss quadrature weights and kronron quadrature abscissae and weights
% as evaluated with 80 decimal digit arithmetic by l. w. fullerton,
% bell labs, nov. 1981.
%
      wg=[ 0.129484966168869693270611432679082e0;
       0.279705391489276667901467771423780e0 ;
       0.381830050505118944950369775488975e0 ;
       0.417959183673469387755102040816327e0 ];
%
      xgk=[ 0.991455371120812639206854697526329e0 ;
       0.949107912342758524526189684047851e0 ;
      0.864864423359769072789712788640926e0 ;
      0.741531185599394439863864773280788e0 ;
      0.586087235467691130294144838258730e0 ;
      0.405845151377397166906606412076961e0 ;
      0.207784955007898467600689403773245e0 ;
      0.000000000000000000000000000000000e0 ];
       
      slope=0.5*(b-a); point=0.5*(b+a); sum=0.0;
      for k=1:4
       ugk=slope*xgk(2*k)+point;
       fval=eval([f,'(ugk)']);
       sum=sum+wg(k)*fval;
      end;
      for k=1:3
      ugk=-slope*xgk(2*k)+point;
       fval=eval([f,'(ugk)']);
       sum=sum+wg(k)*fval;
      end;
      sum=sum*slope;