function plotInternalForces2dTruss(coordinates, connectivity, E, areas, disp, theTitle)
% Create a figure
figure;
hold on;
set(gca, 'Color', 'white'); % Set the axis background color to white
set(gcf, 'Color', 'white'); % Set the figure background color to white

% Plot undeformed structure
for i = 1:size(connectivity, 1)
    % Extract the x and y coordinates of the undeformed structure
    x = [coordinates(connectivity(i, 1), 1), coordinates(connectivity(i, 2), 1)];
    y = [coordinates(connectivity(i, 1), 2), coordinates(connectivity(i, 2), 2)];
    plot(x, y, 'k--', 'Marker', '.','MarkerSize',15);
end

% Compute axial load
N = zeros(size(connectivity, 1), 1);
for i = 1:size(connectivity, 1)
    % Calculate the original length of the element
    l_undeformed = norm(coordinates(connectivity(i,2), :) - coordinates(connectivity(i,1), :));
    theta_undeformed = atan2(coordinates(connectivity(i,2), 2) - coordinates(connectivity(i,1), 2), ...
        coordinates(connectivity(i,2), 1) - coordinates(connectivity(i,1), 1));
    
    % Global reference coordinates
    u_X1 = disp(2 * connectivity(i, 1) - 1);
    u_Y1 = disp(2 * connectivity(i, 1));
    u_X2 = disp(2 * connectivity(i, 2) - 1);
    u_Y2 = disp(2 * connectivity(i, 2));
    
    % Local reference coordinates
    u_x1 = u_X1 * cos(theta_undeformed) + u_Y1 * sin(theta_undeformed);
    u_x2 = u_X2 * cos(theta_undeformed) + u_Y2 * sin(theta_undeformed);
    
    deltaL = u_x2 - u_x1;
    N(i) = E * areas(i) / l_undeformed * deltaL / 1000;
end

% Determine the colormap and normalization
max_abs_N = max(abs(N));
theCMAP='jet';
cmap = colormap(theCMAP);

% Plot internal forces with heat map effect
for i = 1:size(connectivity, 1)
    x = [coordinates(connectivity(i, 1), 1), coordinates(connectivity(i, 2), 1)];
    y = [coordinates(connectivity(i, 1), 2), coordinates(connectivity(i, 2), 2)];
    
    % Normalizing the force for color mapping
    color_idx = round(((N(i) + max_abs_N) / (2 * max_abs_N)) * (size(cmap, 1) - 1)) + 1;
    color_idx = max(1, min(color_idx, size(cmap, 1)));  % Clamp to valid indices
    color = cmap(color_idx, :);
    
    % Plot the element with the corresponding color
    plot(x, y, 'Color', color, 'LineWidth', 2);
    
    % Add the value of N on each element
    mid_x = mean(x);
    mid_y = mean(y);
    text(mid_x, mid_y, sprintf('%.2f', N(i)), 'Color', 'black', 'FontSize', 12, 'HorizontalAlignment', 'center', 'VerticalAlignment', 'middle');
end

% Set axis properties
axis equal;
axis off;

% Add a colorbar
caxis([-max_abs_N, max_abs_N]);
colorbar;
title('Axial Force (kN)');

hold off;
end