% plot Gaussian pdf:
dx=0.1
x = [-4.0:dx:4.0];
p = exp(-0.5*x.^2)/sqrt(2.0*pi);
hold on
plot(x,p, 'linewidth', 10 )

% compute integral between x1 and x2:
x1=1.0
x2=2.0
P = sum(p((x>=x1)&(x<x2)))*dx

% draw random numbers:
r = randn( 10000, 1 );
hist(r,x,1.0/dx)

% check P:
Pr = sum((r>=x1)&(r<x2))/length(r)

hold off