Script has no errors, but no plot is given

3 views (last 30 days)
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
Aquatris
Aquatris on 10 Jul 2024
Edited: Aquatris on 10 Jul 2024
You are not calling the function which you called exp() that produces the plot so why do you expect a plot?
Also never use built-in function/variable names as custom function names. exp() already exist in matlab.
Sergio
Sergio on 10 Jul 2024
Edited: Sergio on 10 Jul 2024
The problem lies in "A(i, i) = -2 / (dr^3) + 3 / (r(i)^2 * dr) + e^(r(i)^2) / (i * hbar^3 / (sqrt(2*m)^3));" I defined the function exp(r ) since it should be connected to the line here. Should I remove the function calling at the top, and rewrite exp(r(i)^2) instead? Did that, then I got a new error related to eigvecs which is unknown.

Sign in to comment.

Accepted Answer

KSSV
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');

More Answers (1)

Mathieu NOE
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');

Categories

Find more on General Physics in Help Center and File Exchange

Tags

Products


Release

R2024a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!