function pt = block_phase_transition_2(N, M)
    K_min = 1;
    K_max = 25;
    K_range = K_min:K_max;
    pt = zeros(25, 1);
    d = N;
    m = M;
    
    for K_idx = 1:length(K_range)
        K = K_range(K_idx);
        s = K;
        syms t;
        syms u;
        f = s*(m+t^2)+(d-s)*int((u-t)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t,inf);
        g = diff(f,t);
        t1 = solve(g);
        n = s*(m+t1^2)+(d-s)*int((u-t1)^2*u^(m-1)*exp(-u^2/2)/(2^(m/2-1)*gamma(m/2)),u,t1,inf);
        pt(K_idx) = n;
    end
end