% ---------------------------
% principal.m
% shows dangers in principal value
% integral computation using 
% evenly spaced points
% ---------------------------
clear all; close all;
N=100001;
x=linspace(-1,1,N);
x0 = 0.1;
h=x(2)-x(1);
% Simpson rule weights (1-3-1 repeated)
w=ones(1,N);
w(2:2:N-1)=4;w(3:2:N-2)=2;
w=w/3;
w = w*h;
int = 0;
for i = 1:N
    if(x(i)~=x0) 
        int = int + w(i).*1./(x(i)-x0);
    else
        disp('singularity skipped')
    end
end
fprintf('%10s %15.5f %10s %15g \n',' x0 = ',x0,' int = ',int);
