[f,t,bcoef,n,k,id,x1,x2,tol,quad,ierr,work]=dbfqad(f,t,bcoef,n,k,id,x1,x2,tol,quad,ierr,work);

function [f,t,bcoef,n,k,id,x1,x2,tol,quad,ierr,work]=dbfqad(f,t,bcoef,n,k,id,x1,x2,tol,quad,ierr,work);
persistent a aa ansmlv b bb iflg il1 il2 ilo inbv left mflag np1 npk q ta tb wtol ; 

if isempty(inbv), inbv=0; end;
%***BEGIN PROLOGUE  DBFQAD
%***PURPOSE  Compute the integral of a product of a function and a
%            derivative of a K-th order B-spline.
%***LIBRARY   SLATEC
%***CATEGORY  H2A2A1, E3, K6
%***TYPE      doubleprecision (BFQAD-S, DBFQAD-D)
%***KEYWORDS  INTEGRAL OF B-SPLINE, QUADRATURE
%***AUTHOR  Amos, D. E., (SNLA)
%***DESCRIPTION
%
%     Abstract    **** a doubleprecision routine ****
%
%         DBFQAD computes the integral on (X1,X2) of a product of a
%         function F and the ID-th derivative of a K-th order B-spline,
%         using the B-representation (T,BCOEF,N,K).  (X1,X2) must be a
%         subinterval of T(K) .LE. X .LE. T(N+1).  An integration rou-
%         tine, DBSGQ8 (a modification of GAUS8), integrates the product
%         on subintervals of (X1,X2) formed by included (distinct) knots
%
%         The maximum number of significant digits obtainable in
%         DBSQAD is the smaller of 18 and the number of digits
%         carried in doubleprecision arithmetic.
%
%     Description of Arguments
%         Input      F,T,BCOEF,X1,X2,TOL are doubleprecision
%           F      - external function of one argument for the
%                    integrand BF(X)=F(X)*DBVALU(T,BCOEF,N,K,ID,X,INBV,
%                    WORK)
%           T      - knot array of length N+K
%           BCOEF  - coefficient array of length N
%           N      - length of coefficient array
%           K      - order of B-spline, K .GE. 1
%           ID     - order of the spline derivative, 0 .LE. ID .LE. K-1
%                    ID=0 gives the spline function
%           X1,X2  - end points of quadrature interval in
%                    T(K) .LE. X .LE. T(N+1)
%           TOL    - desired accuracy for the quadrature, suggest
%                    10.*DTOL .LT. TOL .LE. .1 where DTOL is the maximum
%                    of 1.0D-18 and doubleprecision unit roundoff for
%                    the machine = D1MACH(4)
%
%         Output     QUAD,WORK are doubleprecision
%           QUAD   - integral of BF(X) on (X1,X2)
%           IERR   - a status code
%                    IERR=1  normal return
%                         2  some quadrature on (X1,X2) does not meet
%                            the requested tolerance.
%           WORK   - work vector of length 3*K
%
%     Error Conditions
%         Improper input is a fatal error
%         Some quadrature fails to meet the requested tolerance
%
%***REFERENCES  D. E. Amos, Quadrature subroutines for splines and
%                 B-splines, Report SAND79-1825, Sandia Laboratories,
%                 December 1979.
%***ROUTINES CALLED  D1MACH, DBSGQ8, DINTRV, XERMSG
%***REVISION HISTORY  (YYMMDD)
%   800901  DATE WRITTEN
%   890531  Changed all specific intrinsics to generic.  (WRB)
%   890531  REVISION DATE from Version 3.2
%   891214  Prologue converted to Version 4.0 format.  (BAB)
%   900315  CALLs to XERROR changed to CALLs to XERMSG.  (THJ)
%   900326  Removed duplicate information from DESCRIPTION section.
%           (WRB)
%   920501  Reformatted the REFERENCES section.  (WRB)
%***end PROLOGUE  DBFQAD
%
%
if isempty(iflg), iflg=0; end;
if isempty(ilo), ilo=0; end;
if isempty(il1), il1=0; end;
if isempty(il2), il2=0; end;
if isempty(left), left=0; end;
if isempty(mflag), mflag=0; end;
if isempty(npk), npk=0; end;
if isempty(np1), np1=0; end;
if isempty(a), a=0; end;
if isempty(aa), aa=0; end;
if isempty(ansmlv), ansmlv=0; end;
if isempty(b), b=0; end;
if isempty(bb), bb=0; end;
if isempty(q), q=0; end;
if isempty(ta), ta=0; end;
if isempty(tb), tb=0; end;
if isempty(wtol), wtol=0; end;
t_shape=size(t);t=reshape(t,1,[]);
bcoef_shape=size(bcoef);bcoef=reshape(bcoef,1,[]);
work_shape=size(work);work=reshape(work,1,[]);
%***FIRST EXECUTABLE STATEMENT  DBFQAD
ierr = 1;
quad = 0.0d0;
if( k<1 )
xermsg('SLATEC','DBFQAD','K DOES NOT SATISFY K.GE.1',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
elseif( n<k ) ;
xermsg('SLATEC','DBFQAD','N DOES NOT SATISFY N.GE.K',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
elseif( id<0 || id>=k ) ;
xermsg('SLATEC','DBFQAD','ID DOES NOT SATISFY 0.LE.ID.LT.K',2,1);
else;
[wtol ]=d1mach(4);
wtol = max(wtol,1.0d-18);
if( tol<wtol || tol>0.1d0 )
xermsg('SLATEC','DBFQAD','TOL IS LESS DTOL OR GREATER THAN 0.1',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
else;
aa = min(x1,x2);
bb = max(x1,x2);
if( aa>=t(k) )
np1 = fix(n + 1);
if( bb<=t(np1) )
if( aa==bb )
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
npk = fix(n + k);
%
ilo = 1;
[t,npk,aa,ilo,il1,mflag]=dintrv(t,npk,aa,ilo,il1,mflag);
[t,npk,bb,ilo,il2,mflag]=dintrv(t,npk,bb,ilo,il2,mflag);
if( il2>=np1 )
il2 = fix(n);
end;
inbv = 1;
q = 0.0d0;
for left = il1 : il2;
ta = t(left);
tb = t(left+1);
if( ta~=tb )
a = max(aa,ta);
b = min(bb,tb);
[f,t,bcoef,n,k,id,a,b,inbv,tol,ansmlv,iflg,work]=dbsgq8(f,t,bcoef,n,k,id,a,b,inbv,tol,ansmlv,iflg,work);
if( iflg>1 )
ierr = 2;
end;
q = q + ansmlv;
end;
end; left = fix(il2+1);
if( x1>x2 )
q = -q;
end;
quad = q;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end;
%
%
xermsg('SLATEC','DBFQAD','X1 OR X2 OR BOTH DO NOT SATISFY T(K).LE.X.LE.T(N+1)',2,1);
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
return;
end;
end;
t_shape=zeros(t_shape);t_shape(:)=t(1:numel(t_shape));t=t_shape;
bcoef_shape=zeros(bcoef_shape);bcoef_shape(:)=bcoef(1:numel(bcoef_shape));bcoef=bcoef_shape;
work_shape=zeros(work_shape);work_shape(:)=work(1:numel(work_shape));work=work_shape;
end
%DECK DBHIN