It appears that the default settings for integral2 don't "realize" that the critical region that dominates the integral reisdes in relatively small area far from the origin. This particular problem can be solved by either changing the integration method or doing some work up front to narrow down the integration limits to only cover the region of non-trivial probability density.
Updated code to illustrate both options:
Sigma_bpara = [54.1622939975184 -497.142586464551; -497.142586464551 8986.94991970858];
Mu_bpara = [-0.423193291368079; -934.725277662277];
inputs = [Sigma_bpara(1,1), Sigma_bpara(2,2), Sigma_bpara(1,2), Mu_bpara(1), Mu_bpara(2)];
P1 = Pc_integral(inputs,R,1);
P2 = Pc_integral(inputs,R,2);
function Pc = Pc_integral(inputs,R,option)
Pc_function = @(x,y) (1/(2*pi*sqrt(det(Sigma)))) * exp(-0.5*transpose([x;y]-Mu)*inv(Sigma)*([x;y]-Mu));
Pc=integral2(@(theta,r) arrayfun(pcoll,theta,r), 0,2*pi,0,R,'Method','iterated');
Pc=integral2(@(theta,r) arrayfun(pcoll,theta,r), 3*pi/2-atan(.1),3*pi/2+atan(.1),500,1500);
For option 2 I hard-wired in the limits of integration based a plot of the pdf. In reality, one would have to come up with a procedure to pick those limits based on the mean and covariance.
I don't know enough about the integral2() methods to have an informed opinion on whether or 'iterated' is alwasy satisfactory. My preference is to be more insightful in choosing the limits of integration.