P=20;
$Assumptions=p0≥0&&p1≥0&&p2≥0&&p3≥0&&p0≤P&&p1≤P&&p2≤P&&p3≤P

m0=0;m1=0;m2=0.511;m3=0.511;
en0:=(Sqrt[(p0^2)+(m0^2)]);en1:=(Sqrt[(p1^2)+(m1^2)]);en2:=(Sqrt[(p2^2)+(m2^2)]);(*en3:=(Sqrt[(p3^2)+(m3^2)]);*)
en3:=en0+en1-en2;
p3:=(Sqrt[(en3^2)-(m3^2)]);
region[p0_,p1_,p2_]:=(
en0≥m0&&en1≥m1&&en2≥m2&&en3≥m3&&
Max[p0,p1]≤ Min[p0,p1]+p2+p3&&
Max[p2,p3]≤ Min[p2,p3]+p0+p1
)
RegionPlot3D[region[p0,p1,p2],
{p0,0,P},{p1,0,P},{p2,0,P},
AxesLabel->{"p0","p1","p2"},PlotPoints->20,Mesh->None,
PlotStyle->Directive[Yellow,Opacity[0.5]]]

Manipulate[RegionPlot[
region[p0,p1,p2]/.p1->(Subscript[p, 0]),
{p0,0,P},{p2,0,P}
],{(Subscript[p, 0]),0,P}]

Simplify[region[p0,p1,p2]]

Reduce[Simplify[region[p0,p1,p2]&&$Assumptions]/.{p0->5,p1->7},p2]

NSolve[{
$Assumptions,
p0≥p1&&p2≥p3,
en3≥m3,
p0≤p1+p2+p3,
p2≤p3+p0+p1
},{p2},Reals]

A = 100;a=(10^(-5)); r=1.25;
Solve[A==a((1-(r^n))/(1-r)),n]