function plotElasticModulus(l, nel, sigma, tleng, E0)
    % Generate covariance matrix R using exponential covariance function
    i = 1:nel;
    j = 1:nel;
    ee = zeros(nel, nel);
    R = zeros(nel, nel);

    e(i) = tleng / 2 * nel + (tleng / nel) * (i - 1);
    for i = 1:nel
        for j = 1:nel
            ee(i, j) = e(i) - e(j);
            R(i, j) = sigma^2 * exp(-abs(ee(i, j) / (tleng / nel)));  % Use tleng / nel instead of d1
        end
    end

    % Cholesky decomposition of R
    C1 = chol(R, 'lower');

    % Generate random vectors using Cholesky decomposition
    rng(1, 'twister');
    Z1 = randn(nel, l);
    alpha = C1 * Z1;

    % Elastic modulus calculation
    E = E0 * (1 + alpha);

    % Calculate mean values
    mean1 = mean(E, 2);

    % Plotting
    n0 = tleng / nel;
    n00 = linspace(0, tleng, nel);

    figure;
    plot(n00, E)

    hold on
    line([0, tleng], [E0, E0], 'Color', 'k', 'LineWidth', 2);  % Add a thick black horizontal line at E = E0

    % Identify and mark points where elastic modulus goes below 0
    [r, c] = find(E < 0);
    below_zero_indices = sub2ind(size(E), r, c);  % Convert row and column indices to linear indices
    plot(n00(r), E(below_zero_indices), 'k*', 'MarkerSize', 10, 'LineWidth', 1.5);  % Use 'k*' for a dark star, adjust 'MarkerSize', and add 'LineWidth'
    ylabel('Elastic Modulus in GPa', 'fontweight', 'bold', 'fontsize', 16)  % Include the units for thermal conductivity
    xlabel('Distance from fixed end', 'fontweight', 'bold', 'fontsize', 16)
    set(gca, 'FontSize', 16, 'fontweight', 'bold');
    title('Elastic Modulus Simulation');
end