% wachkyPlot00.m 
% Effects of cancellation and rounding in computing 
% f(x) = (1-cos(x))/x^2
% when x is getting close to sqrt(eps) or less.
% K. Ming Leung, 01/22/04 (scalar vs. vectorized versions)

clear all; 

nPts = 50000; xmax = 1e-7;

tic;
Xvector = linspace(eps,xmax,nPts); % create absissa points
% compute function values in a straightforward (vectorized) way:
Yvector = (1-cos(Xvector))./(Xvector.*Xvector);    
toc;

tic;
del = (xmax-eps)/(nPts-1);
Xvector1(1) = eps;
Yvector1(1) = (1-cos(Xvector1(1)))./(Xvector1(1)^2);    
for n=2:(nPts)
    Xvector1(n) = Xvector(n-1) + del;
    Yvector1(n) = (1-cos(Xvector1(n)))./(Xvector1(n)^2);    
end
toc;

maxAbsDiffX = max(abs(Xvector1-Xvector)),
maxAbsDiffY = max(abs(Yvector1-Yvector)),

plot(Xvector,Yvector,'b-',Xvector1,Yvector1,'r--');
