Script has no errors, but no plot is given

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

[Aquatris]

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]

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.

Answer 1

[KSSV]

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

Answer 2

[Mathieu NOE]

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