@Markiyan Hirnyk 

All the source code is in my post; it is enough a copy & paste + execute.

I had the impression that you have tried lately to adopt a minimal politeness.
It seems that I was wrong. I am sorry.

Edit. It seems also that there is still hope, because you have removed your reply.

@Markiyan Hirnyk 

a = mu*phi-1
b = (1-mu)*phi-1
(The factor not depending on y was removed).

@krismalo 
I don't think it could be simpler; the procedure has only a few lines.
If f has multiple roots, solve(f,x) will find all of them. If you want the minimal one, simply replace
a := solve(f,x);
    with
a := min(solve(f,x) );
 

@Serwa 

In your example there are two sets of functions (obtained by pdsolve), given by "generators":
S1:  f(x,y,z,t) = F1(x, y, z)
S2:  f(x,y,z,t) = F2(x+y, x+z, t)

and you want the "generator" for S1 intersect S2.
You have obtained this generator via pdsolve (for the system of the two PDEs)
but it is not difficult to find it by a simple inspection as:
f(x,y,z,t) = F3(x+y, x+z).
Of course, this generator F3 could appear in another form e.g.
f(x,y,z,t) = F4(x+y, z-y).

If should be possible to write a procedure to find the "generator" for such an intersection
and also to decide whether  2 sets (e.g. F3, F4 above) are equal.
The type of sets given by pdsolve will be needed to write such a procedure.
 

Thank you, Robert.
That's the difference between brute force and a simple but real algorithm (or between O(n^4) and O(n^2)).
Probably if a compiler was not available, many users would search for better algorithms :-)

Best regards,
V. A.

@Robert Israel 

@Markiyan Hirnyk 

Too late for me, and for you too :-)

Needs a serious revision.

Probability(X*Y-Z <= 5);  #???
                               1
Probability(X*Y-Z*U <= 5);
      Error, (in Statistics:-Probability) summand is singular in the interval of summation
Probability(X*Y=0);  #???
                               1/199
Probability(X*Y<1);  #ok
                             19999/39601
Probability(X*Y>-1);
      Error, (in Statistics:-Probability) summand is singular in the interval of summation
Probability(X*Y<0);  #ok
                             19602/39601
Probability(X*Y>0);
      Error, (in Statistics:-Probability) summand is singular in the interval of summation

@Markiyan Hirnyk 

N:=199^4:
pr:=proc()
local a::integer[8],b::integer[8],c::integer[8],d::integer[8],Z::integer[8]:=0;
for a from -99 to 99 do
for b from -99 to 99 do
for c from -99 to 99 do
for d from -99 to 99 do
if a*d-b*c=0 then Z:=Z+1 fi od:od:od:od:
Z;
end:

cpr:=Compiler:-Compile(pr):

cpr():
%/N=evalf(%/N);
561441/1568239201 = 0.3580072476e-3
 

The optimal solution is when |x| + |y| is minimal, where x,y are integers and (17/2)*x+11*y = 6. So:

Optimization:-LPSolve( p+q, {x <=p, x>=-p, y<=q, y>=-q, (17/2)*x+11*y = 6}, assume=integer);
    [3, [p = 2, q = 1, x = 2, y = -1]]

Hence: x = 2, y = -1

@Carl Love 

I am more interested in the following aspect of the ILP.
It is known that unlike simplex, LPSolve uses floats, so the solution is not necessarily exact.
In the case of an ILP with integer coefficients, are there situations/conditions (e.g. on dimensions and magnitude of the coefficients) when the solution can be guaranteed to be exact? Otherwise it would be interesting an example.
Anyway, the integer and binary options for simplex would be very useful.

(BTW, it is known that an ILP is "harder" than a LP wrt complexity.)

@Carl Love 

But wasn't it already solved with "using simplex" below?

Edit. OK, simplex has not an integer option.

@hitstudent 

You should experiment with a familiar surface, in order to understand how it works.

Take e.g. the upper half of a sphere. The following commands produce the same plot:

plot3d(sqrt(1-z^2), theta=0..2*Pi, z=0..1, coords=cylindrical, scaling=constrained);

plot3d( [sqrt(1-z^2), theta, z], theta=0..2*Pi, z=0..1, coords=cylindrical, scaling=constrained);

plot3d( [sqrt(1-z^2)*cos(theta), sqrt(1-z^2)*sin(theta), z], theta=0..2*Pi, z=0..1, scaling=constrained);

@tomleslie 

"Maple 2015 introduced a change in the handling of Pi. Whenever Pi is adjacent to a floating-point number in a sum, product, or power, Pi is automatically converted to a floating-point number.  See Compatibility Issues in Maple 2015.
In order to keep the symbol Pi in such cases, set the kernelopts option floatPi to false."

@taro 

Yes, you have at least a vector plot, but not a professional one.
Here is the plot obtained with the free program Maxima. Compare also the sizes (27 KB vs 1.8 MB).

test1.pdf