Julia: NLsolve with function inside functions

I am trying to solve a system of equations given by the function "func_bcd!". The issue i run into is that Julia "freezes" or keeps running when trying to run the code provided below.

I think the issue is related to how the first equation in the function "func_bcd!" calls two other functions "r_p_this" and "r_m_this". Maybe i am defining the functions in a wrong manner, but i am not sure.

All help much appreciated.

Best regards, Rasmus Damgaard

using Distributions, Random, Cubature, NLsolve, LinearAlgebra;

Random.seed!(123); # Setting the seed

function F(x)
    (1-exp(-x/2))
end;

function G(x)
    min(1,x/C)
end;

function rₘ_this(x)

    function rₘ_this_func(z)
        min(1,x[1]*z[2]./(x[1]*z[1]+x[2]*x[3]).*pdf(Gamma(μᵤ/2,1),z[1]).*pdf(Gamma(μᵤ/2,1),z[2]))
    end;

    (rₘ,rₘ_error) = hcubature(rₘ_this_func,a,b,abstol=1e-8)
    return rₘ
end;

function rₚ_this(x)

    function rₚ_this_func(z)
        min(1,(x[1]*z[2]+x[2]*x[3])./(x[1]*z[1]).*pdf(Gamma(μᵤ/2,1),z[1]).*pdf(Gamma(μᵤ/2,1),z[2]))
    end;

    (rₚ,rₚ_error) = hcubature(rₚ_this_func,a,b,abstol=1e-8)
    return rₚ
end;

function func_bcd!(f,x)
    f[1] = G(1)-G((rₚ_this(x)-rₘ_this(x))/(rₚ_this(x)+rₘ_this(x)))-x[1]
    f[2] = 1-(1-x[2])*x[3]/((1-x[2])*x[3]+(1-G(1))*μᵤ)-rₘ_this(x)
    f[3] = μ_bar*F((1-G(1)*μᵤ*σₚₙ)/(1-x[2])*x[1]+(1-G(1))*μᵤ)-x[1]
end;

a = [0 0];
b = [3*μᵤ 3*μᵤ];
C = 2;
μ_bar = 20;
μᵤ = 60;
σₚₙ = 0.381971;
x₀ = [0.5 0.0 2.92629];

x = nlsolve(func_bcd!, x₀)