% plot Gaussian pdf:
dx=0.01;
x = [-4.0:dx:4.0];
p = exp(-0.5*x.^2)/sqrt(2.0*pi);
plot(x, p, 'linewidth', 4)
hold on
% show area of integral:
x1=1.0;
x2=2.0;
area(x((x>=x1)&(x<=x2)), p((x>=x1)&(x<=x2)), 'FaceColor', 'r' )
hold off

% compute integral between x1 and x2:
P = sum(p((x>=x1)&(x<x2)))*dx;
fprintf('Integral between %.2g and %.2g: %.3g\n', x1, x2, P);

% draw random numbers:
r = randn(10000, 1);

% check P:
Pr = sum((r>=x1)&(r<x2))/length(r);
fprintf('Probability of a number between %.2g and %.2g: %.3g\n', x1, x2, Pr);

% infinite integral:
P = sum(p)*dx;
fprintf('Integral between -infinity and +infinity: %.3g\n', P);