% DDFEX44 Function used in Example 3.44.
%
function y=ddfex44(x1,x2)
lt=find(abs(x1-x2)<1); gt=find(abs(x1-x2)>=1);
if isempty(lt)
  y=(cos(x1)-cos(x2))./(x1-x2);
else
  sx=size(x1);
  y=zeros(sx);
  slt=size(lt,2); y0=zeros(1,slt);
  x=x1(lt)-x2(lt);
  for is=1:slt
    t=x(is); t2=t^2;
    term=.5; s0=0; s=term;
    k=0;
    while s~=s0
      s0=s; k=k+1;
      term=-t2*term/(16*k*(k+.5));
      s=s0+term;
    end
    y0(is)=s;
  end
  y(lt)=-2*sin(.5*(x1(lt)+x2(lt))).*y0;
  y(gt)=(cos(x1(gt))-cos(x2(gt)))./(x1(gt)-x2(gt));
end
