can you add parfor at my code so it run faster:
% MOEMA optimization with enhanced crowding distance and Levy flight
clc; 
clear; 
close all;

% Set problem parameters
NPop = 150;  % Population size
T = 300;     % Number of iterations
dim = 30;    % Number of decision variables
lb = zeros(1, dim);  % Lower bound
ub = ones(1, dim);   % Upper bound

% Run improved MOEMA
[paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, @zdt1);

% Plot the Pareto front
figure;
plot(paretoFront(:,1), paretoFront(:,2), 'ro');
xlabel('f1');
ylabel('f2');
title('Pareto Front for ZDT1 using Improved MOEMA');

% Calculate and display metrics
disp(['Number of solutions in Pareto Front: ', num2str(size(paretoFront, 1))]);
disp(['Minimum f1: ', num2str(min(paretoFront(:,1)))]);
disp(['Minimum f2: ', num2str(min(paretoFront(:,2)))]);

% Plot the convergence of the Pareto front
figure;
for i = 1:T
    front = squeeze(fitness_history(:,i,:));
    [~, unique_idx] = unique(front, 'rows');
    front = front(unique_idx,:);
    plot(front(:,1), front(:,2), 'bo', 'MarkerSize', 2);
    xlabel('f1');
    ylabel('f2');
    title(['Pareto Front Evolution at Iteration ', num2str(i)]);
    xlim([0 1]);
    ylim([0 1]);
    drawnow;
    pause(0.05);
end

% Plot the trajectory of a random solution
figure;
plot(1:T, Trajectories(1,:));
xlabel('Iteration');
ylabel('x1 value');
title('Trajectory of x1 for a Random Non-dominated Solution');

function [paretoFront, paretoSet, Trajectories, fitness_history, position_history] = MOEMA(NPop, T, lb, ub, dim, FITNESSFCN)
    % Preallocation
    fitness_history = zeros(NPop, T, 2);
    position_history = zeros(NPop, T, dim);
    Trajectories = zeros(dim, T);
    
    % Generation of initial population
    pop = initialization(NPop, dim, ub, lb);
    pop_dad = pop(1:NPop/2, :);
    pop_mum = pop(NPop/2+1:end, :);
    
    % Evaluation
    obj_dad = evaluate_population(pop_dad, FITNESSFCN);
    obj_mum = evaluate_population(pop_mum, FITNESSFCN);
    obj = [obj_dad; obj_mum];
    
    % Initialize archive with non-dominated solutions and crowding distances
    archive = initialize_archive(pop, obj);
    
    for gen = 1:T
        k1 = randperm(NPop/2);
        k2 = randperm(NPop/2);
        k1x = pop_dad(k1, :);
        k2x = pop_mum(k2, :);
        
        for k = 1:NPop/2
            % Enhanced mating process with Levy flight
            Imates = 1 + (var(k1x(k, :)) - var(k2x(k, :)));
            levy_step = levy_flight(dim);  % Generate Levy flight step
            p = levy_step ./ max(abs(levy_step));  % Normalize
            
            if Imates >= 0
                temp1 = p .* k1x(k, :);
                temp2 = (1-p) .* k2x(k, :);
                pop_child = temp1 + temp2;
            else
                temp1 = p .* k2x(k, :);
                temp2 = (1-p) .* k1x(k, :);
                pop_child = temp1 + temp2;
            end
            
            % Impact of best environment using crowding distance
            if ~isempty(archive.pop)
                % Select best solution based on crowding distance
                [~, best_idx] = max(archive.crowding_distance);
                best_pop = archive.pop(best_idx, :);
                K = levy_flight(dim) < 0.8;
                pop_child = best_pop .* (1-K) + pop_child .* K;
            end
            
            % Bounding
            pop_child = max(min(ub, pop_child), lb);
            
            % Evaluate child
            obj_child = FITNESSFCN(pop_child);
            
            % Update population using enhanced selection
            [pop_dad, obj_dad] = update_population_with_crowding(pop_dad, obj_dad, pop_child, obj_child, k1(k));
            [pop_mum, obj_mum] = update_population_with_crowding(pop_mum, obj_mum, pop_child, obj_child, k2(k));
            
            % Exploration process with Levy flight
            if rand() < 0.2 && ~isempty(archive.pop)
                levy_step = levy_flight(dim);
                temp1 = archive.pop(randi(size(archive.pop, 1)), :);
                pop_child = temp1 + levy_step;
                pop_child = max(min(ub, pop_child), lb);
                obj_child = FITNESSFCN(pop_child);
                
                [pop_dad, obj_dad] = update_population_with_crowding(pop_dad, obj_dad, pop_child, obj_child, randi(NPop/2));
                [pop_mum, obj_mum] = update_population_with_crowding(pop_mum, obj_mum, pop_child, obj_child, randi(NPop/2));
            end
        end
        
        % Combine pools and update archive
        pop = [pop_dad; pop_mum];
        obj = [obj_dad; obj_mum];
        archive = update_archive_with_crowding(archive, pop, obj);
        
        % Record history
        fitness_history(:, gen, :) = obj;
        position_history(:, gen, :) = pop;
        
        if ~isempty(archive.pop)
            Trajectories(:, gen) = archive.pop(1, :)';
        end
        
        disp(['Iteration #', num2str(gen), ' completed. Archive size: ', num2str(size(archive.pop, 1))]);
    end
    
    paretoFront = archive.obj;
    paretoSet = archive.pop;
end

function step = levy_flight(dim)
    % Implementation of Levy flight
    beta = 1.5;  % Levy distribution parameter
    sigma = (gamma(1+beta)*sin(pi*beta/2)/(gamma((1+beta)/2)*beta*2^((beta-1)/2)))^(1/beta);
    u = randn(1, dim) * sigma;
    v = randn(1, dim);
    step = u./abs(v).^(1/beta);
end

function archive = initialize_archive(pop, obj)
    % Initialize archive with crowding distance
    [is_dominated, ~] = pareto_dominance(obj);
    non_dominated = ~is_dominated;
    
    archive.pop = pop(non_dominated, :);
    archive.obj = obj(non_dominated, :);
    archive.crowding_distance = calculate_enhanced_crowding_distance(archive.obj);
end

function crowding_dist = calculate_enhanced_crowding_distance(objectives)
    [n_solutions, n_objectives] = size(objectives);
    crowding_dist = zeros(n_solutions, 1);
    
    if n_solutions <= 2
        crowding_dist(:) = Inf;
        return;
    end
    
    % Calculate enhanced crowding distance for each objective
    for i = 1:n_objectives
        [sorted_obj, sort_idx] = sort(objectives(:,i));
        
        % Set boundary points to infinity
        crowding_dist(sort_idx(1)) = Inf;
        crowding_dist(sort_idx(end)) = Inf;
        
        obj_range = sorted_obj(end) - sorted_obj(1);
        if obj_range == 0
            continue;
        end
        
        % Enhanced crowding distance calculation
        for j = 2:n_solutions-1
            prev_diff = sorted_obj(j) - sorted_obj(j-1);
            next_diff = sorted_obj(j+1) - sorted_obj(j);
            
            % Weight based on position in objective space
            position_weight = 1 - abs(2*j/n_solutions - 1);
            
            % Enhanced distance considering both spacing and position
            crowding_dist(sort_idx(j)) = crowding_dist(sort_idx(j)) + ...
                (prev_diff + next_diff) * position_weight / obj_range;
        end
    end
end

function archive = update_archive_with_crowding(archive, pop, obj)
    % Combine current archive with new solutions
    combined_pop = [archive.pop; pop];
    combined_obj = [archive.obj; obj];
    
    % Get non-dominated solutions
    [is_dominated, ~] = pareto_dominance(combined_obj);
    non_dominated = ~is_dominated;
    
    % Update archive
    archive.pop = combined_pop(non_dominated, :);
    archive.obj = combined_obj(non_dominated, :);
    
    % Calculate new crowding distances
    archive.crowding_distance = calculate_enhanced_crowding_distance(archive.obj);
    
    % Limit archive size if needed (keep most diverse solutions)
    max_archive_size = 1000;  % Adjust as needed
    if size(archive.pop, 1) > max_archive_size
        [~, sort_idx] = sort(archive.crowding_distance, 'descend');
        archive.pop = archive.pop(sort_idx(1:max_archive_size), :);
        archive.obj = archive.obj(sort_idx(1:max_archive_size), :);
        archive.crowding_distance = archive.crowding_distance(sort_idx(1:max_archive_size));
    end
end

function [pop, obj] = update_population_with_crowding(pop, obj, child, child_obj, index)
    % Create temporary population including the child
    temp_pop = [pop; child];
    temp_obj = [obj; child_obj];
    
    % Calculate crowding distances
    crowding_dist = calculate_enhanced_crowding_distance(temp_obj);
    
    % Non-dominated sorting
    [is_dominated, ~] = pareto_dominance(temp_obj);
    
    if ~is_dominated(end)  % If child is non-dominated
        if is_dominated(index)  % If current solution is dominated
            pop(index, :) = child;
            obj(index, :) = child_obj;
        else  % Both non-dominated, use crowding distance
            if crowding_dist(end) > crowding_dist(index)
                pop(index, :) = child;
                obj(index, :) = child_obj;
            end
        end
    end
end

% Helper Functions

function pop = initialization(NPop, dim, ub, lb)
    pop = rand(NPop, dim) .* (ub - lb) + lb;
end

function obj = evaluate_population(pop, FITNESSFCN)
    obj = zeros(size(pop, 1), 2);  % Assuming 2 objectives
    for i = 1:size(pop, 1)
        obj(i, :) = FITNESSFCN(pop(i, :));
    end
end

function [is_dominated, domination_count] = pareto_dominance(objectives)
    [n_solutions, n_objectives] = size(objectives);
    is_dominated = false(n_solutions, 1);
    domination_count = zeros(n_solutions, 1);

    for i = 1:n_solutions
        for j = i+1:n_solutions
            if dominates(objectives(i, :), objectives(j, :))
                is_dominated(j) = true;
                domination_count(i) = domination_count(i) + 1;
            elseif dominates(objectives(j, :), objectives(i, :))
                is_dominated(i) = true;
                domination_count(j) = domination_count(j) + 1;
            end
        end
    end
end

function result = dominates(obj1, obj2)
    result = all(obj1 <= obj2) && any(obj1 < obj2);
end

function f = zdt1(x)
    f1 = x(1);
    g = 1 + 9 * sum(x(2:end)) / (length(x) - 1);
    h = 1 - sqrt(f1 / g);
    f2 = g * h;
    f = [f1, f2];
end