function [N_b, N_s] = block_phase_transition(N, M)
    tau_min = 0;
    tau_max = 10;
    tau_interval = 0.05;
    tau_range = tau_min:tau_interval:tau_max;
    
    K_min = 1;
    K_max = 15;
    K_range = K_min:K_max;
    
    N_b = zeros(length(K_range), 1);
    N_s = zeros(length(K_range), 1);
    
    cache_filename_1 = "I_2.mat";
    cache_filename_2 = "I_" + string(2 * M) + ".mat";
    
    if exist(cache_filename_1, "file")
        load(cache_filename_1, "I_2");
    else
        I_2 = zeros(length(tau_range), 0);
    
        for tau_idx = 1:length(tau_range)
            tau = tau_range(tau_idx);
            I_2(tau_idx) = calc_block_integral(tau, 2);
        end
        save(cache_filename_1, "I_2");
    end
    
    if exist(cache_filename_2, "file")
        load(cache_filename_2, "I_2M");
    else
        I_2M = zeros(length(tau_range), 0);
    
        parfor tau_idx = 1:length(tau_range)
            tau = tau_range(tau_idx);
            I_2M(tau_idx) = calc_block_integral(tau, 2 * M);
        end
        save(cache_filename_2, "I_2M");
    end
    
    
    for K_idx = 1:length(K_range)
        K = K_range(K_idx);
        f_set_1 = zeros(length(tau_range), 1);
        f_set_2 = zeros(length(tau_range), 1);
    
        for tau_idx = 1:length(tau_range)
            tau = tau_range(tau_idx);
            f_set_1(tau_idx) = 1/2 * (K * (2 * M + tau ^ 2) + (N - K) * I_2M(tau_idx));
            f_set_2(tau_idx) = M/2 * (K * (2 + tau ^ 2) + (N - K) * I_2(tau_idx));
        end
    
        N_b(K_idx) = min(f_set_1);
        N_s(K_idx) = min(f_set_2);
    end
end