How to implement Gaussian Approximation

I'm having some trouble trying to implement the following function:

To express the inverse function of rho I have used a numerical approximation for the second equality of the piecewise defined function, as it does not have an analytical solution.

I have the impression that something must have been done wrongly, as the result should be more 'positive' (greater than 0) for each iteration.

I will leave here my code:

 l = [0 0.3078 0.27287 0 0 0 0.41933];
 r = [0 0 0 0 0 0.4 0.6];

 sigma = 0.8747;

 mu0 = 2/sigma;

 iterations = 50;

 % Density evolution algorithm depiction for finding the treshold of irregular LDPC codes
 syms x;

 l_idle = zeros(1,length(l));
 r_idle = zeros(1,length(r));

 Q_1 = exp(-0.4527*x^0.86 + 0.0218);
 Q_2 = sqrt(pi/x)*exp(-x/4)*(1-20/(7*x));

 mu = zeros(1,iterations+1);

 for k=2:1:iterations+1
    for i = 1:length(l_idle)
        if ((mu0 + (i-1)*mu(k-1)) < 10)
            l_idle(i) = subs(Q_1,x,(mu0 + (i-1)*mu(k-1)));
        else
            l_idle(i) = subs(Q_2,x,(mu0 + (i-1)*mu(k-1)));
        end
    end
    lambda = l*transpose(l_idle);

    for j = 1:length(r_idle)
        b = 1-(1-lambda)^(j-1);
        g = subs(Q_1,x,b);
        if or(subs(Q_1,x,10) < g, subs(Q_1,x,0) >= g)
            r_idle(j) = (1/0.4527*(0.0218-log(g)))^(1/0.86);
        else
            xx = linspace(0, subs(Q_2,x,10), 10000);
            yy = f(xx);
            r_idle(j) = interp1(yy,xx,b);
        end
    end
    mu(k) = r*transpose(r_idle);
end

If anyone has any question, I will be pleased to answer you.

Thank you in advance, and may you have a nice week!