% FEX37  A function used in Examples 3.36 and 3.37. 
%
function y=fex37(x,eta,theta)
n=size(x,1);
y=zeros(n,1);
for k=1:n
  x1=x(k)-eta;
  if abs(x1)<=1
    r=0;
    term=x1;
    s=0; s1=term;
    while s1~=s
      s=s1;
      r=r+1;
      if r>200, disp('series in fex37 does not converge'), return, end
      term=x1*term/(r+1);
      s1=s+term;
    end
%
% The second statement below applies for "large" theta.
%
    y(k)=x(k)*sqrt(1+.5*theta*x(k))/(exp(-x(k))*s);
%     y(k)=x(k)/(exp(-x(k))*s);
  else
%
% The second statement below applies for "large" theta.
%
    y(k)=x(k)*sqrt(1+.5*theta*x(k))/(exp(-eta)-exp(-x(k)));
%    y(k)=x(k)/(exp(-eta)-exp(-x(k)));
  end 
end
