function [z,Y] = qqplot(X)
N = length(X);

Y = sort( (X - mean(X)) / std(X) );
p = (1:N) / (N+1);
z = sqrt(2)*erfinv(2*p - 1);
plot(z, Y, 'o', z, z, 'k');
title('Q-Q plot against normal');
end