Script has no errors, but no plot is given
3 views (last 30 days)
Show older comments
Hi , I wonder what is wrong here with the script. Nothinig happens when it should plot something... Can someone help? Thanks
% Constants
hbar = 1.0545718e-34; % Planck's constant
m = 9.10938356e-31; % Electron mass
l = 1; % angular momentum quantum number
E_ion = 1; % Ionization energy, example value
%Functions
function e =exp(r)
% Discretization parameters
r_min = 0.01; % Avoid division by zero
r_max = 10;
N = 100; % Number of grid points
r = linspace(r_min, r_max, N);
dr = r(2) - r(1);
% Initialize the R(r) vector
R = zeros(N, 1);
% Set up the finite difference matrix
A = zeros(N, N);
for i = 3:N-2
A(i, i-2) = -1 / (2 * dr^3);
A(i, i-1) = 2 / (dr^3) - 3 / (r(i) * dr^2);
A(i, i) = -2 / (dr^3) + 3 / (r(i)^2 * dr) + e^(r(i)^2) / (i * hbar^3 / (sqrt(2*m)^3));
A(i, i+1) = -2 / (dr^3) + 3 / (r(i) * dr^2);
A(i, i+2) = 1 / (2 * dr^3);
end
% Apply boundary conditions (example: R(0) = 0 and R(r_max) = 0)
A(1,1) = 1;
A(N,N) = 1;
% Solve the eigenvalue problem
[~, D] = eig(A);
% The eigenvalues are the diagonal elements of D
eigenvalues = diag(D);
% The solution R(r) corresponds to the eigenvector with eigenvalue closest to E_ion
[~, idx] = min(abs(eigenvalues - E_ion));
R = eigvecs(:, idx);
% Plot the solution
plot(r, R);
xlabel('r');
ylabel('R(r)');
title('Radial Wavefunction');
end
2 Comments
Accepted Answer
KSSV
on 10 Jul 2024
% Constants
hbar = 1.0545718e-34; % Planck's constant
m = 9.10938356e-31; % Electron mass
l = 1; % angular momentum quantum number
E_ion = 1; % Ionization energy, example value
%Functions
% Discretization parameters
r_min = 0.01; % Avoid division by zero
r_max = 10;
N = 100; % Number of grid points
r = linspace(r_min, r_max, N);
dr = r(2) - r(1);
% Initialize the R(r) vector
R = zeros(N, 1);
% Set up the finite difference matrix
A = zeros(N, N);
for i = 3:N-2
A(i, i-2) = -1 / (2 * dr^3);
A(i, i-1) = 2 / (dr^3) - 3 / (r(i) * dr^2);
A(i, i) = -2 / (dr^3) + 3 / (r(i)^2 * dr) + exp(r(i)^2) / (i * hbar^3 / (sqrt(2*m)^3));
A(i, i+1) = -2 / (dr^3) + 3 / (r(i) * dr^2);
A(i, i+2) = 1 / (2 * dr^3);
end
% Apply boundary conditions (example: R(0) = 0 and R(r_max) = 0)
A(1,1) = 1;
A(N,N) = 1;
% Solve the eigenvalue problem
[eigvecs, D] = eig(A);
% The eigenvalues are the diagonal elements of D
eigenvalues = diag(D);
% The solution R(r) corresponds to the eigenvector with eigenvalue closest to E_ion
[~, idx] = min(abs(eigenvalues - E_ion));
R = eigvecs(:, idx);
% Plot the solution
plot(r, R);
xlabel('r');
ylabel('R(r)');
title('Radial Wavefunction');
0 Comments
More Answers (1)
Mathieu NOE
on 10 Jul 2024
again...
let's fix it
I cannot check here if e^(r(i)^2) (probably wrong) must be replaced by exp(r(i)^2) - I let you double check that point , but that is actually what I did
then, you have to change one more line to define eigvecs
[~, D] = eig(A); => [eigvecs, D] = eig(A);
% Constants
hbar = 1.0545718e-34; % Planck's constant
m = 9.10938356e-31; % Electron mass
l = 1; % angular momentum quantum number
E_ion = 1; % Ionization energy, example value
% Discretization parameters
r_min = 0.01; % Avoid division by zero
r_max = 10;
N = 100; % Number of grid points
r = linspace(r_min, r_max, N);
dr = r(2) - r(1);
% Initialize the R(r) vector
R = zeros(N, 1);
% Set up the finite difference matrix
A = zeros(N, N);
for i = 3:N-2
A(i, i-2) = -1 / (2 * dr^3);
A(i, i-1) = 2 / (dr^3) - 3 / (r(i) * dr^2);
A(i, i) = -2 / (dr^3) + 3 / (r(i)^2 * dr) + exp(r(i)^2) / (i * hbar^3 / (sqrt(2*m)^3));
A(i, i+1) = -2 / (dr^3) + 3 / (r(i) * dr^2);
A(i, i+2) = 1 / (2 * dr^3);
end
% Apply boundary conditions (example: R(0) = 0 and R(r_max) = 0)
A(1,1) = 1;
A(N,N) = 1;
% Solve the eigenvalue problem
[eigvecs, D] = eig(A);
% The eigenvalues are the diagonal elements of D
eigenvalues = diag(D);
% The solution R(r) corresponds to the eigenvector with eigenvalue closest to E_ion
[~, idx] = min(abs(eigenvalues - E_ion));
R = eigvecs(:, idx);
% Plot the solution
plot(r, R);
xlabel('r');
ylabel('R(r)');
title('Radial Wavefunction');
0 Comments
See Also
Categories
Find more on General Physics in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!