[x y] = meshgrid([-10:0.25:10]);

n = 25;

N = 1/(2*pi)^(3/4)*1/sqrt(factorial(2*n));

psi = N^2*(x.^2+y.^2).^n.*exp(-1/2*(x.^2+y.^2));
surf(x,y,psi);
xlabel('x');
ylabel('y');
title('|\Psi|^2');