function z = experfc (x,yy);
% computes product z = exp(x)*erfc(yy) in a manner that does not blow up
% based on rational Chebyshev approximations (Wexler, 1989, USGS OFR 89-56)
p1=zeros(1,5);
q1=zeros(1,5);
p2=zeros(1,9);
q2=zeros(1,9);
p3=zeros(1,6);
q3=zeros(1,6);
p1=[3.20937795E03 3.774852376E02 1.386415415E02...
    3.16112374E00 1.857777062E-01];
q1=[2.84423683E03 1.282616526E03 2.440246379D02...
    2.360129095E01 1.0];
p2=[1.230339355E03 2.051078378E03 1.712047613E03...
   8.819522212E02 2.986351382E02 6.611919064E01...
   8.883149794E00 5.6418849699E-01 2.1531153547E-08];
q2=[1.230339355E03 3.439367674144E03 4.3626190901E03...
   3.2907992357E03 1.621389574567E03 5.3718110186E02...
   1.1769395089E02 1.574492611071E01 1.0];
p3=[-6.567491615E-04 -1.60837851487E-02 -1.25781726111229E-01...
   -3.603448999498E-01 -3.053266349612D-01 -1.631538713730E-02];
q3=[2.3352049762687E-03 6.051834131244E-02 5.279051029514E-01...
    1.8729528499234E00 2.5685201922898E00 1.0];

if yy==0.0;
  z=exp(x);
end;

if yy <= 0.46875;
  sump=0.0;
  sumq=0.0;
  for ii=1:1:5;
     y2i=yy^(2*ii-1);
     sump=sump+p1(ii)*y2i;
     sumq=sumq+q1(ii)*y2i;
   end;
erf=yy*sump/sumq;
        if yy < 0.0;
        erf=-erf;
        end;
erfcy=1.0-erf;
z=exp(x)*erfcy;
end;

if yy > 0.46875;
sump=0.0;
sumq=0.0;
  for ii=1:1:9;
  	yi=yy^(ii-1);
  	sump=sump+p2(ii)*yi;
  	sumq=sumq+q2(ii)*yi;
  end;
z=exp(x-yy*yy)*sump/sumq;
   if yy < 0.0;
   z=2.0*exp(x)-z;
   end;
end;

if yy > 4.0;
sump=0.0;
sumq=0.0;
  for ii=1:1:6;
  y2i=yy^(-2*(ii-2));
  sump=sump+p3(ii)*y2i;
  sumq=sumq+q3(ii)*y2i;
end;
sqrtpi=0.5641895835477562869481;
z=sqrtpi+sump/(yy*yy*sumq);
z=exp(x-yy*yy)*z/yy;
  if yy < 0.0;
  z=2.0*exp(x)-z;
  end;
end;





