% wackyPlotSol.m 
% Effects of cancellation and rounding in computing 
% f(x) = (1-cos(x))/x^2  for x around sqrt(eps) or less.
% By directly computing the locations of the discontinuities,
% we numerically verify that our interpretation is correct.
% K. Ming Leung, 01/25/04

clear all; hold off;
x=linspace(eps,1e-7,50000);
x2 = x.*x;  % x^2
y0 = (1-cos(x))./(x2);  % straightforward function computation
plot(x,y0,'b-','LineWidth',2.5); 
hold on;                % want to add more data to the graph     
title('f(x) = (1 - cos(x))/x^2','fontSize',16);
xlabel('x','fontSize',24);
ylabel('f(x)','fontSize',24);
set(gca,'FontSize',14);

% plot the jumps at the discontinuities
plotNumber = 10;        % number of discontinuities to plot
eMach = eps/2;          % machine epsilon
idx = 0:1:plotNumber;   % 0, 1, 2, ..., plotNumber
odd = 2*idx+1;          % 1, 3, 5, ..., 2*plotNumber+1
xn = sqrt(odd*eMach);   % location of discontinuities
lower = idx./odd;       % f(x) right before a jump
upper = (idx+1)./odd;   % f(x) right after a jump
plot(xn,lower,'or',xn,upper,'or');

x22 = 0.5*x2;
yBest = 0.5*(sin(x22)./x22).^2;
plot(x,yBest,'g-','LineWidth',1.0);