I try to evaluate the following integral in MATLAB: $$\int_0^\infty f_R(R)$$ in which it's expected to be 1. Where we have: $$f_R(R)=\frac d{dR}\left(\exp\left(-2\pi\lambda_m\left(\int_0^Rre^{-\beta r}\ dr+\int_0^{R^{\frac ab}}r(1-e^{-\beta r})\ dr\right)\right)\right)$$ For example by setting $a=2,b=4,\lambda_m = 5\times 10^{-6}, \beta = 0.0143 $. But the result is $236.7063$!!! Here is my code:
a = 2;
b = 4;
Req = @(r) r^(a/b);
bei = 0.0143;
lambdamm = 5e-6;
Plr = @(r) r*exp(-bei *r);
Plr2 = @(r) r*(1-exp(-bei*r));
iner1 = @(R) integral(Plr,0,R,'ArrayValued',true);
iner2 = @(R) integral(Plr2,0,Req(R),'ArrayValued',true);
Lambdar2 = @(r) (exp(-2$\pi$ * lambdamm*iner1(r)))*exp(-2$\pi$ * lambdamm * iner2(r));
lambdar2 = @(r) 2$\pi$ * lambdamm * ((r * Plr(r))+((a/b) * (r^(a/b))) * (1-Plr(r)));
f_Lr = @(r) Lambdar2(r) * lambdar2(r);
integral(f_Lr,0,inf,'ArrayValued',true)
This is not an answer but it is too long for a comment (and moreover I do not know Matlab).
I do not understand why you define in the post $$f_R(R)=\frac{d\phi(R)}{dR}$$ and then need to integrate $f_R(R)$.
To me, $$\int f_R(R)\,dR=\phi(R)$$ and then the problem is to compute $$I=\int_0^\infty f_R(R)\,dR=\lim_{R\to \infty } \, \phi (R)-\lim_{R\to 0 } \, \phi (R)$$
Since all the calculations can be done analytically, I checked using your values, and the result is $-1$ which does not seem to agree with your expectation.