% ----------------------------------------------------
% principal3.m
% Principal value
% integral computation using 
% Riemann sum; points symmetrically wrt. singularity
% not vectorized
% integral sin(x^2)/(x-x0) from -1 to 1
% ----------------------------------------------------
clear all; close all;
h=10^-5;
Nx=2/h;
x0 = -0.95;
res=[];
while x0<1
    int = 0;        
    for i=1:Nx
        x=x0+i*h;
        if(x<1)
            int = int + sin(x.^2)/(x-x0);
        end
        x=x0-i*h;
        if(x>-1) 
            int = int + sin(x.^2)/(x-x0);
        end
    end
    int = int*h;
    res=[res [x0 int]' ];
    fprintf('%10s %15.5f %10s %15.5f \n',' x0 = ',x0,' int = ',int);
    x0 = x0 + 0.05;
end
plot(res(1,:),res(2,:),'bx');
legend('numerical');
xlabel('x0');
