Plotting contours of a probability density

How do you create a script file to plot contours of a probability density in the y=0 plane? This is what I have done, but it's obviously wrong. Please help!

z=[-6:0.01:6];

x=[-6:0.01:6];

r=sqrt(x.^2+z.^2);

u1=(8*pi)^(-0.5)*(1-(r./2)).*exp(-r./2);

u2=(8*pi)^(-0.5)*(r./2).*(z/r).*exp(-r./2);

v1=(u1-u2).^(2).*0.5;

contour(x,z,v1)

xlabel('x')

zlabel('z')

I have attached the plot I am supposed to get.

Answer 1

[Star Strider]

You’re missing a meshgrid call and some element-wise dot-operator division

Otherwise, your code is correct:

z=[-6:0.01:6];
x=[-6:0.01:6];  
[X,Z] = meshgrid(x,z);  
r=sqrt(X.^2+Z.^2);  
u1=sqrt(8*pi)*(1-(r./2)).*exp(-r./2);  
u2=sqrt(8*pi)*(r./2).*(Z./r).*exp(-r./2);  
v1=(u1-u2).^(2).*0.5;

figure(1)
contour(X,Z,v1,20)
xlabel('x')
zlabel('z')

This isn’t exactly like the plot you posted, but it’s close! I’ll leave it to you to supply the necessary refinements to get it looking the way you want. You may want to adjust the number of contours the plot draws (I opted for 20 here).

[Niki]

Thank you! I didn't know about adding the number of contours and I was not sure about adding meshgrid.

[Star Strider]

My pleasure!

Answer 2

[Youssef Khmou]

You have to use two dimensional arrays of x and z to produce two dimensional pdf, try :

   [x,z]=meshgrid(-6:0.01:6);
r=sqrt(x.^2+z.^2);
u1=(8*pi)^(-0.5)*(1-(r./2)).*exp(-r./2);
u2=(8*pi)^(-0.5)*(r./2).*(z./r).*exp(-r./2);
v1=(u1-u2).^(2).*0.5;
contour(x,z,v1)
xlabel('x')
zlabel('z')

[Niki]

Thanks for your prompt reply! Very useful.