function y = tanh(xin)
%TANH   Hyperbolic tangent.
%   TANH(X) is the hyperbolic tangent of the elements of X.

%   C. Moler, 6-17-92
%   Copyright (c) 1984-98 by The MathWorks, Inc.
%   $Revision$  $Date$

%   21 Oct 92 (jmh) -- fixed for matrix inputs
%        2 Sep 93 (jmh) -- fixed for null matrix input

% Complex argument
if ~isreal(xin) %  2 Sep 93 (jmh)
    tr = tanh(real(xin));
    ti = i*tan(imag(xin));
    y = (tr + ti)./(1 + tr.*ti);

% Real argument
else

   % Reference: W. J. Cody and W. Waite, "Software Manual
   % for the Elementary Functions", 1980, chapter 13.
   
   xbig = 18.715;                       % 0.5*log(4/eps)
   xmid = 0.5493061443340549;           % log(3)/2
   xsmall = 1.49e-8;                    % sqrt(eps)
   x = abs(xin);
   y = zeros(size(x));

   % tanh(x) == 1 to working precision
   k = find(xbig <= x);
   if ~isempty(k)           % 21 Oct 92 (jmh)
      y(k) = ones(size(find(k)));       % added FIND -- 20 Jul 92 (jmh)
   end

   % 1/2 < tanh(x) < 1
   k = find((xmid < x) & (x < xbig));
   if ~isempty(k)           % 21 Oct 92 (jmh)
      y(k) = 1 - 2./(exp(2*x(k))+1);
   end

   % x < tanh(x) <= 1/2
   k = find((xsmall < x) & (x <= xmid));
   if ~isempty(k)           % 21 Oct 92 (jmh)
      p = [-0.16134119023996228e4 -0.99225929672236083e2 -0.96437492777225470];
      q = [0.48402357071988689e4 0.22337720718962313e4  0.11274474380534949e3];
      xx = x(k).^2;
      y(k) = x(k) + x(k).*(xx.*(((p(3)*xx + p(2)).*xx + p(1))) ./ ...
                          (((xx + q(3)).*xx + q(2)).*xx + q(1)));
   end

   % tanh(x) == x to working precision
   k = find((x <= xsmall) | isnan(x));
   if ~isempty(k)       % 21 Oct 92 (jmh)
      y(k) = x(k);
   end

   y = sign(xin).*y;
end
