eq11=1- ( (1-Rab)*(1-(Rbc*Rac)))==1;
eq22=1- ( (1-Rbc)*(1-(Rab*Rac)))==1;
eq33=1- ( (1-Rac)*(1-(Rbc*Rab)))==1;
Let us look at those equations for a moment. eq11, eq22, eq33 are independent of x.
All three of them are equations, and you multiply the equations by expressions that involve x. The multiplication acts on both sides of the equation, so if you had f(x) == c and multiply both sides by constant A then A*f(x) == A*c would be the resulting equation.
The multiplication does not introduce new solutions: the resulting extended equation has solutions either when the original equations (eq11, eq22, eq33) are true by themselves, or when the multipliers (gam1, gam2, gam3) are 0. So because eq11, eq22, eq33 are independent of x, the multiplied values, eq1, eq2, eq3, will have solutions at the points where solve([eq11, eq22, eq33], [Rab, Rbc, Rac]) works, and will also have solutions where solve([gam1, gam2, gam3], x) works.
S = solve(eq1,eq2,eq3,Rab,Rbc,Rac);
And there you see the solutions that are independent of x, as discussed above. You could plot them over x, but it would be constant graphs.
Okay, how about the case where the x contributions become 0, solve([gam1, gam2, gam3],x) ?
You can see that for each individual gam* that most of the solutions are imaginary, and all of them are outside x in 0:56000 . You can also see that although the first of them shares a real solution with the third, and the second of them shares a real solution with the third, that there is no real-valued x for which all three have a solution. Therefore when you solve the three of them as a group, you are not going to see any real-valued solutions.
... Not unless solve() has omitted some solution. So let us try:
G = sum([gam1,gam2,gam3].^2);
Wait, is that a solution?
Looks like it??
... But no, it is a false solution, where the values reach 0 because of round-off, but are non-zero. See
Gsol2 = vpasolve(simplify(G));
More digits, the solution extends out a lot further... in fact there is no solution until you reach x = infinity.
function p = Gamcdf(x,a,b)
function y = Gammainc(x,a)
y = int(t^(a-1)*exp(-t), t, 0, x) / gamma(a);