Code covered by the BSD License

### Highlights fromslatec

from slatec by Ben Barrowes
The slatec library converted into matlab functions.

[f,a,b,result,abserr,resabs,resasc]=qk15(f,a,b,result,abserr,resabs,resasc);
```function [f,a,b,result,abserr,resabs,resasc]=qk15(f,a,b,result,abserr,resabs,resasc);
%***BEGIN PROLOGUE  QK15
%***PURPOSE  To compute I = Integral of F over (A,B), with error
%                           estimate
%                       J = integral of ABS(F) over (A,B)
%***CATEGORY  H2A1A2
%***TYPE      SINGLE PRECISION (QK15-S, DQK15-D)
%***AUTHOR  Piessens, Robert
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%           de Doncker, Elise
%             Applied Mathematics and Programming Division
%             K. U. Leuven
%***DESCRIPTION
%
%           Integration rules
%           Standard fortran subroutine
%           Real version
%
%           PARAMETERS
%            ON ENTRY
%              F      - Real
%                       function subprogram defining the integrand
%                       function F(X). The actual name for F needs to be
%                       Declared E X T E R N A L in the calling program.
%
%              A      - Real
%                       Lower limit of integration
%
%              B      - Real
%                       Upper limit of integration
%
%            ON RETURN
%              RESULT - Real
%                       Approximation to the integral I
%                       Result is computed by applying the 15-POINT
%                       KRONROD RULE (RESK) obtained by optimal addition
%                       of abscissae to the 7-POINT GAUSS RULE(RESG).
%
%              ABSERR - Real
%                       Estimate of the modulus of the absolute error,
%                       which should not exceed ABS(I-RESULT)
%
%              RESABS - Real
%                       Approximation to the integral J
%
%              RESASC - Real
%                       Approximation to the integral of ABS(F-I/(B-A))
%                       over (A,B)
%
%***REFERENCES  (NONE)
%***ROUTINES CALLED  R1MACH
%***REVISION HISTORY  (YYMMDD)
%   800101  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)
%***end PROLOGUE  QK15
%
persistent absc centr dhlgth epmach fc firstCall fsum fv1 fv2 fval1 fval2 hlgth j jtw jtwm1 resg resk reskh uflow wg wgk xgk ; if isempty(firstCall),firstCall=1;end;

if isempty(absc), absc=0; end;
if isempty(centr), centr=0; end;
if isempty(dhlgth), dhlgth=0; end;
if isempty(epmach), epmach=0; end;
if isempty(fc), fc=0; end;
if isempty(fsum), fsum=0; end;
if isempty(fval1), fval1=0; end;
if isempty(fval2), fval2=0; end;
if isempty(fv1), fv1=zeros(1,7); end;
if isempty(fv2), fv2=zeros(1,7); end;
if isempty(hlgth), hlgth=0; end;
if isempty(resg), resg=0; end;
if isempty(resk), resk=0; end;
if isempty(reskh), reskh=0; end;
if isempty(uflow), uflow=0; end;
if isempty(wg), wg=zeros(1,4); end;
if isempty(wgk), wgk=zeros(1,8); end;
if isempty(xgk), xgk=zeros(1,8); end;
if isempty(j), j=0; end;
if isempty(jtw), jtw=0; end;
if isempty(jtwm1), jtwm1=0; end;
%
%
%           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
%
if firstCall,   xgk(1) =[0.9914553711208126e+00];  end;
if firstCall,  xgk(2) =[0.9491079123427585e+00];  end;
if firstCall,  xgk(3) =[0.8648644233597691e+00];  end;
if firstCall,  xgk(4) =[0.7415311855993944e+00];  end;
if firstCall,  xgk(5) =[0.5860872354676911e+00];  end;
if firstCall,  xgk(6) =[0.4058451513773972e+00];  end;
if firstCall,  xgk(7)=[0.2077849550078985e+00];  end;
if firstCall,  xgk(8)=[0.0e+00];  end;
if firstCall,   wgk(1) =[0.2293532201052922e-01];  end;
if firstCall,  wgk(2) =[0.6309209262997855e-01];  end;
if firstCall,  wgk(3) =[0.1047900103222502e+00];  end;
if firstCall,  wgk(4) =[0.1406532597155259e+00];  end;
if firstCall,  wgk(5) =[0.1690047266392679e+00];  end;
if firstCall,  wgk(6) =[0.1903505780647854e+00];  end;
if firstCall,  wgk(7)=[0.2044329400752989e+00];  end;
if firstCall,  wgk(8)=[0.2094821410847278e+00];  end;
if firstCall,   wg(1) =[0.1294849661688697e+00];  end;
if firstCall,  wg(2) =[0.2797053914892767e+00];  end;
if firstCall,  wg(3) =[0.3818300505051189e+00];  end;
if firstCall,  wg(4)=[0.4179591836734694e+00];  end;
firstCall=0;
%
%
%           LIST OF MAJOR VARIABLES
%           -----------------------
%
%           CENTR  - MID POINT OF THE INTERVAL
%           HLGTH  - HALF-LENGTH OF THE INTERVAL
%           ABSC   - ABSCISSA
%           FVAL*  - FUNCTION VALUE
%           RESG   - RESULT OF THE 7-POINT GAUSS FORMULA
%           RESK   - RESULT OF THE 15-POINT KRONROD FORMULA
%           RESKH  - APPROXIMATION TO THE MEAN VALUE OF F OVER (A,B),
%                    I.E. TO I/(B-A)
%
%           MACHINE DEPENDENT CONSTANTS
%           ---------------------------
%
%           EPMACH IS THE LARGEST RELATIVE SPACING.
%           UFLOW IS THE SMALLEST POSITIVE MAGNITUDE.
%
%***FIRST EXECUTABLE STATEMENT  QK15
[epmach ]=r1mach(4);
[uflow ]=r1mach(1);
%
centr = 0.5e+00.*(a+b);
hlgth = 0.5e+00.*(b-a);
dhlgth = abs(hlgth);
%
%           COMPUTE THE 15-POINT KRONROD APPROXIMATION TO
%           THE INTEGRAL, AND ESTIMATE THE ABSOLUTE ERROR.
%
fc = f(centr);
resg = fc.*wg(4);
resk = fc.*wgk(8);
resabs = abs(resk);
for j = 1 : 3;
jtw = fix(j.*2);
absc = hlgth.*xgk(jtw);
fval1 = f(centr-absc);
fval2 = f(centr+absc);
fv1(jtw) = fval1;
fv2(jtw) = fval2;
fsum = fval1 + fval2;
resg = resg + wg(j).*fsum;
resk = resk + wgk(jtw).*fsum;
resabs = resabs + wgk(jtw).*(abs(fval1)+abs(fval2));
end; j = fix(3+1);
for j = 1 : 4;
jtwm1 = fix(j.*2 - 1);
absc = hlgth.*xgk(jtwm1);
fval1 = f(centr-absc);
fval2 = f(centr+absc);
fv1(jtwm1) = fval1;
fv2(jtwm1) = fval2;
fsum = fval1 + fval2;
resk = resk + wgk(jtwm1).*fsum;
resabs = resabs + wgk(jtwm1).*(abs(fval1)+abs(fval2));
end; j = fix(4+1);
reskh = resk.*0.5e+00;
resasc = wgk(8).*abs(fc-reskh);
for j = 1 : 7;
resasc = resasc + wgk(j).*(abs(fv1(j)-reskh)+abs(fv2(j)-reskh));
end; j = fix(7+1);
result = resk.*hlgth;
resabs = resabs.*dhlgth;
resasc = resasc.*dhlgth;
abserr = abs((resk-resg).*hlgth);
if( resasc~=0.0e+00 && abserr~=0.0e+00 )
abserr = resasc.*min(0.1e+01,(0.2e+03.*abserr./resasc).^1.5e+00);
end;
if( resabs>uflow./(0.5e+02.*epmach) )
abserr = max((epmach.*0.5e+02).*resabs,abserr);
end;
end
%DECK QK15I
```