I am trying to obtain $\mathsf{P}_{n}$, the probability that $n$ nodes appear in the ``intersection'' of two circles (with the same radius $r$) formed by two separate nodes, whose area is denoted as $\mathsf{A}$. Notice that the position of each node follows an i.i.d. uniform distribution; the area $\mathsf{A}$ is a function of the distance between the two nodes, $\mathsf{l}$; accordingly, $\mathsf{A}$ becomes a random variable. Although I was able to obtain the pdf for $\mathsf{A}$, since it is not found in a closed form, I am trying to numerically sum the integral given below. \begin{align}\label{eq_pn0} \mathsf{P}_{n} &= \mathbb{E}_{\mathsf{A}} \bigg[ \mathbb{P}\big[ \mathbb{N}\left[ \text{# nodes in } \mathsf{A} \right] = n \big] \bigg]\nonumber\\ &\stackrel{(a)}{=} \mathbb{E}_{\mathsf{A}} \left[ e^{ -\lambda \mathsf{A} }\frac{\left( \lambda \mathsf{A} \right)^{n}}{n!} \right]\nonumber\\ &= \frac{\lambda^{n}}{n!} \displaystyle \int_{\min {\mathsf{A}}}^{\max {\mathsf{A}}} e^{ -\lambda \mathsf{a} } \mathsf{a}^{n} f_{\mathsf{A}}\left(\mathsf{a}\right) \text{d}\mathsf{a} \end{align} (a) follows from the fact that placement of nodes follows a Poisson point process (PPP).
Question: The values for the parameters that I used are $\min \mathsf{A} = 10^6 \text{ m}^2$, $\max \mathsf{A} = 3 \times 10^6 \text{ m}^2$, and $\lambda = 0.01 / m^2$ (meaning 1 node per 10 $\text{m}$ $\times$ 10 $\text{m}$).
Here is the problem. With these values, $\lambda \mathsf{A}$ becomes way too much large, yielding $e^{-\lambda \mathsf{A}}$ to be computed as 0 in MATLAB. This, in turn, leads the entire integral to all zero, regardless the value of $n$. (I know it does not technically exactly zero, and only means an even small number than the computer's precision.) But, based on intuition, the probability of having a certain number of nodes in the area $\mathsf{A}$ certainly should not be that much small!
I still think I am using realistic values for $\min \mathsf{A}$, $\max \mathsf{A}$, and $\lambda$. Am I misconceiving in dealing with the ``unit'' between the intensity $\lambda$ and area $\mathsf{A}$? If not, how did other people ever solved this similar problem?

a) Given your problem formulation, the Poisson assumption is not the only valid point process. Just because the random points are "uniformly" located, does NOT mean that they belong to a (homogeneous) Poisson point process.
Generally speaking, if the points are uniformly located and each point is independently located to all other points, then the points usually belong to either 1) a homogeneous Poisson point process OR 2) a binomial point process. More specifically, in the whole simulation/sample window (or, in mathematical terminology, carrier space or state space) W:
1) if the total number of points is a Poisson variable, then it's a (homogeneous) Poisson point process; 2) if the total number of points is a fixed number m, then it's a binomial point process.
(If you condition on the Poisson point process to have n points, then you get a binomial point process.)
Then the number of points in a subregion/subset A of the simulation window W is a Poisson random variable (for case 1) or a binomial random variable (for case 2).
https://en.wikipedia.org/wiki/Poisson_distribution
https://en.wikipedia.org/wiki/Binomial_distribution
To simulate these two point processes on a disk in MATLAB, Python or R, see these two respective posts:
https://hpaulkeeler.com/point-process-simulation/
https://hpaulkeeler.com/poisson-point-process-simulation/
To simulate on a disk, see this post:
https://hpaulkeeler.com/simulating-a-poisson-point-process-on-a-disk/
The code is also found here: https://github.com/hpaulkeeler/posts
b) The answer doesn't involve an integral.
Let's write the area of the whole simulation/sample window as |W|. Similarly, the area of the subregion is |A|. Then the parameters of two above random variables are (I will use the notation from the above Wikipedia articles, but with primes, so if Wikipedia writes n, I'll write n'):
1) lambda'=lambda * |A|. That is to say, the mean/parameter of the Poisson variable is lambda * |A|. 2) p'=|A|/|W|, q'=1- |A|/|W| , n'=m, k'=n. That is to say, you're flipping $m$ times a biased coin, which has a probability |A|/|W| of a head appearing, and then you're asking what is the probability of n heads appearing?
Now you have to insert the numbers into the expressions for the probabilities (ie evaluate the probability mass functions).
Of course, in practice, for basically all numerical cases, the Poisson point process so closely approximates a binomial point process, one always uses a Poisson point process, because it's so much easier to work with, while losing essentially no accuracy.
In the Poisson case, the final answer is:
P_n=(|A|lambda)^n exp(-|A|lambda)/n!
If you want to know the area of two overlapping circles, then you can, for example, use this site:
http://mathworld.wolfram.com/Circle-CircleIntersection.html
It also appears (equation 3.5) in this paper:
"A stochastic analysis of a greedy routing scheme in sensor networks" by Keeler and Taylor.
To support your result with simulations, the easiest thing would be to simulate the point process on a rectangle that bounds (ie covers) the two disks. Then you can use the equation of a disk (in Cartesian coordiates) to check the points are in a disk. That is, for a disk with radius r centered at (a,b), a point (x,y) exists in the disk if the following (MATLAB expression) holds:
For more information, there are many sources on the related fields of spatial statistics and stochastic geometry. I would recommend the (Springer published) lectures "Spatial Point Processes and their Applications" by Adrian Baddeley. They are found in this book. You can read the entire here book with a Springer subscription.
Actually, anything by Baddeley, including various papers and books, is worth looking at. These lectures also cover the two point processes, but they might be a bit too mathematical. Also see the references at the end of the two posts I mentioned before.
I hope this helps.
PS: Here's some MATLAB code for comparing the probabilities under the binomial and Poisson distributions. Of course, you can also write code to simulate the two point processes.