function [E1] = myexpint(z)

  if (nargin != 1)
    print_usage ();
  endif

  sz = size(z); # Extract the size
  z = z(:); # Write the input as column
  E1 = zeros(size(z)); # Initialize the result
  # In the following, we extract the three zones in which we divide the input
  flag_s = (((real(z) + 19.5).^2./(20.5^2) + imag(z).^2./(10^2)) <= 1) | (real(z)<0 & abs(imag(z))<=1e-08);
  flag_cf = ((((real(z) + 1).^2./(38^2) + imag(z).^2./(40^2)) <= 1) & (~flag_s)) & (real(z) <= 35);
  flag_b = ~flag_s & ~flag_cf;
  range_s = find(flag_s);
  range_cf = find(flag_cf);
  range_b = find(flag_b);
  z_s = z(range_s);
  z_cf = z(range_cf);
  z_b = z(range_b);
  # Series expansion, Abramowitz Stegun
  gm =  0.577215664901532860606512090082402431042159335;
  E1_s = -gm - log(z_s);
  res = -z_s;
  ssum = res;
  k = 1;
  fflag = ones(size(res));
  rrange = find(fflag);
  while(any(fflag))
    rrange = rrange(find(fflag));
    res_tmp = res(rrange);
    res_tmp .*= k*(-z_s(rrange))/((k + 1)^2);
    ssum(rrange) += res_tmp;
    k++;
    res(rrange) = res_tmp;
    fflag = abs(res_tmp) > eps*abs(ssum(rrange));
  endwhile
  E1_s -= ssum;
  # Continued fraction, modified Lentz’s algorithm
  f_new = 2^(-100)*ones(size(z_cf));
  C_new = f_new;
  D_new = zeros(size(z_cf));
  k = 0;
  fflag = ones(size(z_cf));
  rrange = find(fflag);
  Delta = ones(size(z_cf));
  C_old = ones(size(z_cf));
  D_old = ones(size(z_cf));
  f_old = ones(size(z_cf));
  while (k < 1e3 && any(fflag))
    rrange = rrange(find(fflag));
    f_new_tmp = f_new(rrange);
    C_new_tmp = C_new(rrange);
    D_new_tmp = D_new(rrange);
    f_old = f_new_tmp;
    C_old = C_new_tmp;
    D_old = D_new_tmp;
    b = z_cf(rrange)*(mod(k,2) == 0) + (mod(k,2) == 1);
    a = exp(-z_cf(rrange))*(k == 0) + ceil((k)/2)*(k >= 1);
    D_new_tmp = b + a.*D_old;
    D_new_tmp = D_new_tmp.*(D_new_tmp ~= 0) + 2^(-100)*(D_new_tmp == 0);
    C_new_tmp = b + a./C_old;
    C_new_tmp = C_new_tmp.*(C_new_tmp ~= 0) + 2^(-100)*(C_new_tmp == 0);
    D_new_tmp = 1./D_new_tmp;
    Delta(rrange) = C_new_tmp.*D_new_tmp;
    f_new(rrange) = f_old.*Delta(rrange);
    C_new(rrange) = C_new_tmp;
    D_new(rrange) = D_new_tmp;
    fflag = abs(Delta(rrange)-1) > 2^(-50);
    k++;
  endwhile
  E1_cf = f_new;
  # Asymptotic series, Bleistein Handelsman
  Inf_range = find(real(z_b) == -Inf);
  z_b(Inf_range) = -realmax + imag(z_b(Inf_range)); # corrector when real(z)=-Inf
  E1_b = exp(-z_b)./z_b;
  res = ones(size(z_b));
  oldres = res;
  ssum = res;
  k = 0;
  fflag = ones(size(z_b));
  rrange = find(fflag);
  while (k < 1e3 && any(fflag))
    rrange = rrange(find(fflag));
    res_tmp = res(rrange);
    oldres_tmp = oldres(rrange);
    oldres_tmp = res_tmp;
    res_tmp ./= (-z_b(rrange)/(k+1));
    ssum(rrange) += res_tmp;
    k++;
    res(rrange) = res_tmp;
    oldres(rrange) = oldres_tmp;
    fflag = abs(oldres_tmp) > abs(res_tmp);
  endwhile
  E1_b .*= ssum;
  # Re-arranging
  E1(range_s) = E1_s;
  E1(range_cf) = E1_cf;
  E1(range_b) = E1_b;
  E1 = reshape(E1,sz);
endfunction

## Test
# The following values were computed with the octave symbolic package
%!test
%! X =[-50 - 50i  -30 - 50i  -10 - 50i    5 - 50i   15 - 50i   25 - 50i
%!  -50 - 30i  -30 - 30i  -10 - 30i    5 - 30i   15 - 30i   25 - 30i
%!  -50 - 10i  -30 - 10i  -10 - 10i    5 - 10i   15 - 10i   25 - 10i
%!  -50 +  5i  -30 +  5i  -10 +  5i    5 +  5i   15 +  5i   25 +  5i
%!  -50 + 15i  -30 + 15i  -10 + 15i    5 + 15i   15 + 15i   25 + 15i
%!  -50 + 25i  -30 + 25i  -10 + 25i    5 + 25i   15 + 25i   25 + 25i];
%! y_exp = [ -3.61285286166493e+19 + 6.46488018613387e+19i  -4.74939752018180e+10 + 1.78647798300364e+11i   3.78788822381261e+01 + 4.31742823558278e+02i 5.02062497548626e-05 + 1.23967883532795e-04i   3.16785290137650e-09 + 4.88866651583182e-09i   1.66999261039533e-13 + 1.81161508735941e-13i
%!   3.47121527628275e+19 + 8.33104448629260e+19i   1.54596484273693e+11 + 2.04179357837414e+11i   6.33946547999647e+02 + 3.02965459323125e+02i 2.19834747595065e-04 - 9.25266900230165e-06i   8.49515487435091e-09 - 2.95133588338825e-09i   2.96635342439717e-13 - 1.85401806861382e-13i
%!   9.65535916388246e+19 + 3.78654062133933e+19i   3.38477774418380e+11 + 8.37063899960569e+10i   1.57615042657685e+03 - 4.33777639047543e+02i 2.36176542789578e-05 - 5.75861972980636e-04i  -6.83624588479039e-09 - 1.47230889442175e-08i  -2.93020801760942e-13 - 4.03912221595793e-13i
%!  -1.94572937469407e+19 - 1.03494929263031e+20i  -4.22385087573180e+10 - 3.61103191095041e+11i   4.89771220858552e+02 - 2.09175729060712e+03i 7.26650666035639e-04 + 4.71027801635222e-04i   1.02146578536128e-08 + 1.51813977370467e-08i   2.41628751621686e-13 + 4.66309048729523e-13i
%!   5.42351559144068e+19 + 8.54503231614651e+19i   1.22886461074544e+11 + 3.03555953589323e+11i  -2.13050339387819e+02 + 1.23853666784218e+03i -3.68087391884738e-04 + 1.94003994408861e-04i  -1.39355838231763e-08 + 6.57189276453356e-10i  -4.55133112151501e-13 - 8.46035902535333e-14i
%!  -7.75482228205081e+19 - 5.36017490438329e+19i  -1.85284579257329e+11 - 2.08761110392897e+11i  -1.74210199269860e+02 - 8.09467914953486e+02i 9.40470496160143e-05 - 2.44265223110736e-04i   6.64487526601190e-09 - 7.87242868014498e-09i   3.10273337426175e-13 - 2.28030229776792e-13i];
%! assert (y_exp,myexpint(X),-1e-12)

## Exceptional values (-Inf, Inf, NaN, 0, 0.37250741078)
%!test
%! x = [-Inf; Inf; NaN; 0; -0.3725074107813668];
%! y_exp = [-Inf - i*pi; 0; NaN; Inf;0 - i*pi];
%! y = myexpint (x);
%! assert (y, y_exp, 5*eps);

## Test input validation
%!error myexpint ()
%!error myexpint (1,2)
