function optimizeHeatExchanger101()
    % NSGA-II parameters
    nPop = 1000; % Population size
    maxGen = 30; % Maximum number of generations
    pCrossover = 0.999; % Crossover probability

    % Optimization options
    options = gaoptimset('PopulationType', 'doubleVector', ...
                         'PopulationSize', nPop, ...
                         'Generations', maxGen, ...
                         'CrossoverFraction', pCrossover, ...
                         'ParetoFraction', 0.35, ...
                         'PlotFcn', @gaplotpareto);

    % Lower and upper bounds for design variables
    lb = [5, 3, 1, 0.15, 300];
    ub = [10, 6, 2, 0.3, 800];

    % Run the NSGA-II optimizer
    [x, fval] = gamultiobj(@plateFinHeatExchanger, 5, [], [], [], [], lb, ub, options);

    % Display the results to the Command Window
    disp('Optimized Design Variables:');
    disp(x);
    disp('Objective Function Values:');
    disp(fval);
    
    % Obtain the design parameters at which the optimum values of the objective functions are achieved
    designParametersAtOptimum = findOptimalDesignParameters(x, fval);
    disp('Design Parameters at Optimal Objective Function Values:');
    disp(designParametersAtOptimum);
    
    % Plot the objective function values vs Reynolds number
    plotObjectiveFunctionsVsReynolds(x, fval);
    
    % Plot the j/f ratio vs Reynolds number
    plotJFRatioVsReynolds(x, fval);
    
    % Plot the objective function values vs design parameters h, l, s, and t
    plotObjectiveFunctionsVsDesignParameters(x, fval);
end

function plotObjectiveFunctionsVsDesignParameters(x, fval)
    designParams = x(:, 1:4);
    f1_values = fval(:, 1);
    f2_values = fval(:, 2);

    figure;
    
    subplot(2, 2, 1);
    plot(designParams(:, 1), f1_values, 'b.');
    xlabel('h');
    ylabel('Objective Function 1 Value (Negative Colburn Factor)');
    title('Objective Function 1 vs h');

    subplot(2, 2, 2);
    plot(designParams(:, 2), f1_values, 'r.');
    xlabel('l');
    ylabel('Objective Function 1 Value (Negative Colburn Factor)');
    title('Objective Function 1 vs l');
    
    subplot(2, 2, 3);
    plot(designParams(:, 3), f1_values, 'g.');
    xlabel('s');
    ylabel('Objective Function 1 Value (Negative Colburn Factor)');
    title('Objective Function 1 vs s');
    
    subplot(2, 2, 4);
    plot(designParams(:, 4), f1_values, 'm.');
    xlabel('t');
    ylabel('Objective Function 1 Value (Negative Colburn Factor)');
    title('Objective Function 1 vs t');
    
    figure;
    
    subplot(2, 2, 1);
    plot(designParams(:, 1), f2_values, 'b.');
    xlabel('h');
    ylabel('Objective Function 2 Value (Friction Factor)');
    title('Objective Function 2 vs h');
    
    subplot(2, 2, 2);
    plot(designParams(:, 2), f2_values, 'r.');
    xlabel('l');
    ylabel('Objective Function 2 Value (Friction Factor)');
    title('Objective Function 2 vs l');
    
    subplot(2, 2, 3);
    plot(designParams(:, 3), f2_values, 'g.');
    xlabel('s');
    ylabel('Objective Function 2 Value (Friction Factor)');
    title('Objective Function 2 vs s');
    
    subplot(2, 2, 4);
    plot(designParams(:, 4), f2_values, 'm.');
    xlabel('t');
    ylabel('Objective Function 2 Value (Friction Factor)');
    title('Objective Function 2 vs t');
end

% (Remaining functions remain unchanged)
function plotObjectiveFunctionsVsReynolds(x, fval)
    ReynoldsNumbers = x(:, 5);
    f1_values = fval(:, 1);
    f2_values = fval(:, 2);

    figure;
    
    subplot(2, 1, 1);
    plot(ReynoldsNumbers, f1_values, 'b.');
    xlabel('Reynolds Number');
    ylabel('Objective Function 1 Value (Negative Colburn Factor)');
    title('Objective Function 1 vs Reynolds Number');

    subplot(2, 1, 2);
    plot(ReynoldsNumbers, f2_values, 'r.');
    xlabel('Reynolds Number');
    ylabel('Objective Function 2 Value (Friction Factor)');
    title('Objective Function 2 vs Reynolds Number');
end

function plotJFRatioVsReynolds(x, fval)
    ReynoldsNumbers = x(:, 5);
    jf_ratio = fval(:, 1) ./ fval(:, 2);

    figure;
    plot(ReynoldsNumbers, jf_ratio, 'g.');
    xlabel('Reynolds Number');
    ylabel('j/f Ratio');
    title('j/f Ratio vs Reynolds Number');
end

function [f, c] = plateFinHeatExchanger(x)
    h = x(1);
    l = x(2);
    s = x(3);
    t = x(4);
    Re = x(5);

    f1 = -ColburnFactor(h, l, s, t, Re);
    f2 = FrictionFactor(h, l, s, t, Re);
    
    f = [f1, f2]; 
    c = [];
end

function j = ColburnFactor(h, l, s, t, Re)
    if (Re >= 300 && Re <= 800)
        j = 0.661 * (Re^(-0.651)) * ((s/h)^(-0.343)) * ((t/s)^(-0.538)) * ((t/l)^(0.305));
    elseif (Re >= 1000 && Re <= 15000)
        j = 0.185 * (Re^(-0.396)) * ((s/h)^(-0.178)) * ((t/s)^(-0.403)) * ((t/l)^(0.29));
    else
        error('Reynolds number is out of the valid range for the provided formula.');
    end
end

function f = FrictionFactor(h, l, s, t, Re)
    if (Re >= 300 && Re <= 800)
        f = 10.882 * (Re^(-0.79)) * ((s/h)^(-0.359)) * ((t/s)^(-0.187)) * ((t/l)^(0.284));
    elseif (Re >= 1000 && Re <= 15000)
        f = 2.237 * (Re^(-0.236)) * ((s/h)^(-0.347)) * ((t/s)^(-0.151)) * ((t/l)^(0.639));
    else
        error('Reynolds number is out of the valid range for the provided formula.');
    end
end

function designParameters = findOptimalDesignParameters(x, fval)
    [~, minIdx] = min(fval, [], 1);
    designParameters = x(minIdx, :);
end