function [M,dM] = KummerMz(a,b,z);
% 02.06.2014

%  This function computes the Kummer confluent hypergeometric function M(a,b,z)
%  by using the method described by Abad and Sesma 

%%  'z' can be a vector


%%  ATTENTION!!!   Problem for 2*b-4*a = 0


% dM is the derivative dM/da

M = zeros(size(z));
dM = zeros(size(z));


J = find(abs(z*a) < 1);
if not(isempty(J))    % standard computation 
  tol = 1e-6;
  n = 1;  nmax = 50;  mm = length(z(J));
  xn = a*z(J)/b;  M(J) = ones(1,mm);
  while (n < nmax) && (max(abs(xn)) > tol)
    M(J) = M(J) + xn; 
    xn = xn .* z(J)*(a+n)/(b+n)/(n+1);
    n = n+1;  
  end
end

J = find(abs(z*a) >= 1);
if not(isempty(J))
  x = sqrt(z(J)*(2*b-4*a));  mm = length(x);

  n = 1;  nmax = 10;
  fk = zeros(nmax+1);   gk = zeros(nmax+1,mm);  Bk = zeros(nmax+1);
  Pn = zeros(nmax+1,mm);   Pn(1,:) = ones(1,mm);   
  Qn = zeros(nmax+1,mm);   Qn(2,:) = ones(1,mm);
  fk(1) = 1;  gk(1,:) = ones(1,mm);

  M1 = ones(1,mm);  M2 = zeros(1,mm); 
  dM1 = zeros(1,mm);  dM2 = ones(1,mm); 
  while (n < nmax) 
    Bk(n) = bernoulli(2*n);    
    for j=0:n-1,
      fk(n+1) = fk(n+1) + binomial(2*n-1,2*j)*4^(n-j)*abs(Bk(n-j)) * fk(j+1)/(n-j);
    end
    fk(n+1) = fk(n+1)*(1-b/2);
    
    for j=0:floor((n-1)/2),  
      gk(n+1,:) = gk(n+1,:) + binomial(n-1,2*j)*4^(j+1)*abs(Bk(j+1)) * gk(n-2*j,:)/(j+1);
    end
    gk(n+1,:) = gk(n+1,:) .* (-(1i)*z(J)/4);
    
    pn = zeros(1,mm);
    for k=0:floor(n/2),
      pn = pn + binomial(n,2*k) * fk(k+1) * gk(n-2*k+1,:);
    end
    pn = pn .* (1i*z(J)).^n/factorial(n);

    Pn(n+2,:) = 2*(b-1+n)*Pn(n+1,:)./x - Pn(n,:);
    Qn(n+2,:) = 2*(b-1+n)*Qn(n+1,:)./x - Qn(n,:);

    xn1 = pn .* Pn(n+1,:)./x.^n;  xn2 = pn .* Qn(n+1,:)./x.^(n-1);
    dxn1 = pn .* Pn(n+2,:)./x.^(n+1);  dxn2 = pn .* Qn(n+2,:)./x.^n;
    M1 = M1 + xn1;   M2 = M2 + xn2;
    dM1 = dM1 + dxn1;   dM2 = dM2 + dxn2;
    n = n + 1;
  end
  M(J) = gamma(b)* 2^(b-1) * exp(z(J)/2) .* (M1.* x.^(1-b) .* BesselJ(b-1,x) + M2 .* x.^(-b) .* BesselJ(b,x));
  dM(J) = gamma(b)* 2^(b) * z(J).*exp(z(J)/2) .* (dM1.* x.^(1-b) .* BesselJ(b-1,x) + dM2 .* x.^(-b) .* BesselJ(b,x));
end
  

  