Untitled

% Constants
a = 3e-10;         % Lattice constant in meters
b = 3e-10;         % Potential width in meters
V0 = 1.60218e-19;  % Potential height in Joules
hbar = 1.0545718e-34;  % Reduced Planck's constant in J·s
m = 9.10938356e-31;    % Electron mass in kg

% Define a "safe" sine function that avoids division by zero.
% For very small arguments (|x| < tol), sin(x)/x is approximated by 1.
tol = 1e-8;
safe_sin = @(x) (abs(x) < tol).*1 + (abs(x) >= tol).*(sin(x)./x);

% Define the range of k values (crystal momentum)
k = linspace(-pi/a, pi/a, 1000);

% Initialize arrays to store energies for two bands
E1 = NaN(1, length(k));  % First band energy
E2 = NaN(1, length(k));  % Second band energy

% Loop over each k value
for i = 1:length(k)
    % For the current k, define the function f(k1) that must be zero:
    % f(k1) = cos(k(i)*a) - (m*V0*b*a/hbar^2)*[safe_sin(k1*a)] - cos(k1*a)
    func = @(k1) cos(k(i)*a) - (m*V0*b*a/hbar^2)*safe_sin(k1*a) - cos(k1*a);
    
    % Define a range for k1 (avoid k1 = 0 to steer clear of division by zero)
    k1_range = linspace(0.1, 2*pi/a, 100);  % 100 points from 0.1 to 2*pi/a
    
    k1_solutions = [];  % To store found roots for this k
    
    % Bracket the root by checking consecutive intervals
    for j = 1:length(k1_range)-1
        x1 = k1_range(j);
        x2 = k1_range(j+1);
        f1 = func(x1);
        f2 = func(x2);
        % Only proceed if both endpoint values are finite...
        if ~isfinite(f1) || ~isfinite(f2)
            continue;
        end
        % ...and there is a sign change (or one of the endpoints is zero)
        if f1 * f2 <= 0  
            try
                root = fzero(func, [x1, x2]);
                % Add the root if it is finite and not already in the list (within tolerance)
                if isfinite(root)
                    if isempty(k1_solutions) || ~any(abs(k1_solutions - root) < 1e-8)
                        k1_solutions = [k1_solutions, root];
                    end
                end
            catch
                continue;
            end
        end
    end
    
    % Sort the roots and, if at least two are found, record the two lowest-energy solutions.
    k1_solutions = sort(k1_solutions);
    if length(k1_solutions) >= 2
        E1(i) = (hbar^2 * k1_solutions(1)^2) / (2*m); % First band energy
        E2(i) = (hbar^2 * k1_solutions(2)^2) / (2*m); % Second band energy
    end
end

% Calculate overall band statistics using only valid k values (where both bands were found)
validIndices = ~isnan(E1) & ~isnan(E2);
k_valid = k(validIndices);
E1_valid = E1(validIndices);
E2_valid = E2(validIndices);

% Find the minimum and maximum energies for the first two bands (in Joules)
E1_min = min(E1_valid);
E1_max = max(E1_valid);
E2_min = min(E2_valid);
E2_max = max(E2_valid);

% Calculate band widths (in Joules)
band_width_1 = E1_max - E1_min;
band_width_2 = E2_max - E2_min;

% Calculate energy gap between the first and second bands
energy_gap = max(0, E2_min - E1_max);

% Display results in electron-volts (eV)
fprintf('Width of the first energy band: %.4f eV\n', band_width_1 / 1.60218e-19);
fprintf('Width of the second energy band: %.4f eV\n', band_width_2 / 1.60218e-19);
fprintf('Energy gap between the first and second bands: %.4f eV\n', energy_gap / 1.60218e-19);

% -------------------------
% Plot the energy bands vs. k
% Convert energies from Joules to eV for plotting.
E1_eV = E1_valid / 1.60218e-19;
E2_eV = E2_valid / 1.60218e-19;

figure;
plot(k_valid, E1_eV, 'b-', 'LineWidth', 2);
hold on;
plot(k_valid, E2_eV, 'r-', 'LineWidth', 2);
xlabel('Crystal Momentum k (1/m)');
ylabel('Energy (eV)');
title('Energy Bands');
legend('Band 1', 'Band 2', 'Location', 'Best');
grid on;