I wrote a code - it only works for a simple equation. For a more complex one below, it doesn't!
I'm trying to observe how the solution for p3 changes with small 'additions' to the equation.
Code working on simple equation
tabletry = ( Table[p1*p2, {p1, 0, 3}, {p2, 0, 5}] ) // MatrixForm
tabletest = ( Map[ {# == 0} &, tabletry, {-1} ] ) // MatrixForm
deltap = 10^-2;
eqn1[p1_?NumericQ, p2_?NumericQ, p3_] := 1/p3 - 5 p3 == p1 + p2
eqn2[p1_?NumericQ, p2_?NumericQ, p3_] := 1/p3 - 5 p3 + deltap == p1 + p2
diff[p1_?NumericQ,
p2_?NumericQ] := (p3 /. FindRoot[eqn2[p1, p2, p3], {p3, 1}]) - (p3 /. FindRoot[eqn1[p1, p2, p3], {p3, 1}])
(Array[# #2 /. {(0) -> Quiet@Style[Evaluate[diff[#, #2]], Red], _ :> 0} &, {4, 6}, 0]) // MatrixForm
Code NOT working on more complex equation
Variables
C1 = 10^(-1);
C2 = 0.1*C1;
R = 50;
Tb = 0.1;
Geb = 5.;
Z0 = 50;
L[Te_] := 1. + 1.*(Te - 0.1);
Zlcr[Te_, w_] := (1/R + 1/(I*L[Te]*w) + I*C1*w)^-1;
Zload[Te_, w_] := -I*w*C2 + Zlcr[Te, w];
\[CapitalGamma][Te_, w_] := (Zload[Te, w] - Z0)/(Zload[Te, w] + Z0);
y[Te_, w_] := (Abs[\[CapitalGamma][Te, w]])^2;
p[Te_, w_] := Abs[\[CapitalGamma][Te, w]]
DeltaPlocal = 10.^-5;
Table of values
tt2 = Table[ Te /. Chop@Solve[ (1 - y[Te, w]) Pprobe == (Te - Tb) Geb, Te , Method -> Reduce] ,
{w, 2.5, 3., 0.1} , {Pprobe, 1., 2.5, 0.1} ];
tttt2 = ( Map[Chop@(Max[#] - Min[#] &), tt2, {2}] ) // MatrixForm
testmatrix = Map[ {# == 0} &, tttt2, {-1} ]
eqn8[w_?NumericQ, Pprobe_?NumericQ, Te_] := (1 - y[Te, w]) Pprobe == (Te - Tb) Geb
eqn9[w_?NumericQ, Pprobe_?NumericQ, Te_] := DeltaPlocal + (1 - y[Te, w]) Pprobe == (Te - Tb) Geb
diff2[w_?NumericQ, Pprobe_?NumericQ] := ((Te /. FindRoot[eqn9[w, Pprobe, Te], {Te, 0.2}]) - (Te /. FindRoot[eqn8[w, Pprobe, Te], {Te, 0.2}]))
(Array[# #2 /. {(0) -> Quiet@Style[Evaluate[diff2[#, #2]], Red], _ :> 0} &, {6, 16}, 0]) // MatrixForm