% Hypergeometric functions used in Table 4-1 and Table 4-2
function [iF, F, dF, FF, xspan] = Fsph(A,x)
  x=x(:);
  xLim=0.99;
  j=find(x <= xLim);
  if ~isempty(j)
    Sum1=0;
    dSum1=0;
    iSum1=0;
    M=1e7;
    G=1;
    for m=1:M
      g=(A-m+m^2)/m^2;
      G=G*g;
      iS1=1/(m+1)*x(j).^m*G;
      S1=x(j).^m*G;
      dS1=m*x(j).^(m-1)*G;
      iSum1=iSum1+iS1;
      Sum1=Sum1+S1;
      dSum1=dSum1+dS1;
      if max(max(S1/Sum1))<1e-10 && max(max(iS1/iSum1))<1e-10
        break;
      end
    end
    if m == M  warning = 'series did not converge'
    end
    iF(j)=1+iSum1;
    F(j)=1+Sum1;
    dF(j)=dSum1;
  end
  j=find(x > xLim);
  if ~isempty(j)
    [iF0 F0 dF0] = Fsph(A,xLim);
    xspan=fliplr(1-logspace(-9,log10(.01),501));
    xspan(end)=1;
    dFF=@(F,x)[F(2); F(3); (A*F(2)-(1-2*x)*F(3))/x/(1-x)];
    [FF, xspan] = RK4(dFF, [xLim*iF0 F0 dF0], xspan);
    iF(j) = interp1(xspan,FF(1,:),x(j))./x(j);
    F(j) = interp1(xspan,FF(2,:),x(j));
    dF(j) = interp1(xspan,FF(3,:),x(j));
  end
end
