It looks like you're encountering an error because surf expects X, Y, and Z to be 2D arrays, but you are passing 3D arrays. Additionally, the way you are computing Psi seems to be incorrect as you are using a 1D loop over a 3D grid.
Check this code that should fix the issue you encountered:
% Constants
hbar = 1.0545718e-34;
m = 9.10938356e-31;
E_ion = 5.139 * 1.60218e-19;
k_f = 2 * m * E_ion / hbar^2;
% Define alpha (renamed to avoid conflict with MATLAB function)
alpha_val = sqrt(k_f);
% Radial wave function
function R = radial_wavefunction(r, n, l, alpha)
L = laguerreL(n-l-1, 2*l+1, alpha * r.^2);
R = sqrt((2 * alpha)^(l+1) / factorial(n-l-1)) .* exp(-alpha * r.^2) .* (alpha * r).^l .* L;
end
% Spherical harmonic (assuming it's defined elsewhere)
function Y = spherical_harmonic(theta, phi, l, m)
Y = legendre(l, cos(theta)) .* exp(1i * m * phi);
end
% Total wave function in spherical coordinates
function psi = spherical_wavefunction(r, theta, phi, n, l, m, alpha)
R = radial_wavefunction(r, n, l, alpha);
Y = spherical_harmonic(theta, phi, l, m);
psi = R .* Y;
end
% Define grid
r = linspace(0, 10, 50); % Radial coordinate r
theta = linspace(0, pi, 50); % Polar angle theta
phi = linspace(0, 2*pi, 50); % Azimuthal angle phi
% Create grid for 3D plotting
[R, Theta, Phi] = ndgrid(r, theta, phi);
n = 1;
l = 0;
m = 0;
% Compute the wave function on the grid
Psi = abs(spherical_wavefunction(R, Theta, Phi, n, l, m, alpha_val)).^2;
% Spherical to Cartesian Conversion
X = R .* sin(Theta) .* cos(Phi);
Y = R .* sin(Theta) .* sin(Phi);
Z = R .* cos(Theta);
% Plotting 3D surface
figure;
surf(X(:,:,1), Y(:,:,1), Z(:,:,1), Psi(:,:,1)); % Plotting a slice of the 3D data
xlabel('x');
ylabel('y');
zlabel('z');
title(['|\psi_{', num2str(n), ',', num2str(l), ',', num2str(m), '}(r, \theta, \phi)|^2 for Sodium']);
colorbar;
axis equal;
Changed meshgrid to ndgrid to properly handle 3D grids. Used ndgrid to generate a 3D grid and compute the wave function on this grid. Since surf expects 2D inputs, I plotted a slice of the 3D data (Psi(:,:,1)).
I hope this helps!
