function y = betacore(x, a, b)
%BETACORE Core algorithm for the incomplete beta function.
%   BETACORE(x,a,b) is used twice by BETAINC(X,A,B).
%
%   See also BETAINC.

%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision: 5.4 $  $Date: 1997/11/21 23:45:14 $

y = x;
qab = a + b;
qap = a + 1;
qam = a - 1;
am = ones(size(x));
bm = am;
y = am;
bz = 1 - qab .* x ./ qap;
d = zeros(size(x));
app = d;
ap = d;
bpp = d;
bp = d;
yold = d;
m = 1;
while any(abs(y(:)-yold(:)) > 10*eps*abs(y(:)))
   tem = 2 * m;
   d = m * (b - m) .* x ./ ((qam + tem) .* (a + tem));
   ap = y + d .* am;
   bp = bz + d .* bm;
   d = -(a + m) .* (qab + m) .* x ./ ((a + tem) .* (qap + tem));
   app = ap + d .* y;
   bpp = bp + d .* bz;
   yold = y;
   am = ap ./ bpp;
   bm = bp ./ bpp;
   y = app ./ bpp;
   if m == 1
      bz = 1;   % only need to do this first time through
   end
   m = m + 1;
end