It might be that the OP could get by with discerning whether the expression depends on the suffixed variables found.

Download suff_dep.mw

And if one does not prefer type `satisfies` then the type `suffixed` result could be sieved for `depends` by using `select`.

@NIMA112 In your latest worksheet it seems that phi[1] = 2*phi[2] along x=0. That would explain why the contours would not match there, since here phi[1] <> phi[2].

You could adjust by a factor of two, as in this attachment.
h1_h2_ok_ac_2020_3_ac.mw

Also, sometimes you might get a contour mismatch when the close surfaces are quite flat near a specific contour. In that case you could omit the problem contour (or possibly? increase working precision Digits enough to get the needed accuracy) to avoid the appparent mistmatch. In one of your examples in this sheet it seems to suffice to use an even number of contours, thus avoiding having a contour in the (problematic) very middle of the range.

@sursumCorda That difference (in whether the second argument to solve is a set) is a very useful controller. There have been previous Questions involving simplification of relations in which it was the key element. 

And with the result being a set of relations, one could be very slightly terser here:

     `and`(solve(r1 or r2,{x})[])

For fun, one could also handle the equivalent RealRange's using OrProp, instead of throwing it all to solve.

In general I might first turn the active and into inert And, to reduce the risk of inadvertant evaluation to true/false, etc. 

Download ineq_OrProp.mw

@NIMA112 

h1_h2_ok_ac_2020.2.mw

Did it work as described, for anyone interested?

I have corrected the "Product", from Maple 2018 to Maple 18 in which your attachment was last saved.

Those two versions are different, released in different years.

(No need to respond. I know that you are using Maple 18.)

@C_R The evalf inside F is not actually required there. The commands plot, Maximize, and fsolve will evalf themselves.
slow_maximize_and_fsolve_ac_noevalf.mw

You also asked why the following did not work (without substitution y=Y, or :-y=y if you use lowercase y as the name of the parameter of F).

F := y->evalf(Int(unapply(op(1,ig),x),rhs(op(2,ig))));

The y that is the parameter of procedure F is not the same as the global name y already present in op(1,ig).

Even the first of these next two examples has two different "y" names in play.

The command unapply exists to make it easier to build procedures from given expressions.

In the variants under discussion, I was building a procedure that returns an Int call containing a procedure (as integrand). And the original form of the integrand is a given existing expression.

Building procedures that return procedures (ie. a "procedure factory"), substituting for different variables at each level, is pretty awkward without unapply.

The first pair of these look understandable. But that last one, lacking all unapply, is more klunky. Without any unapply, things can be even klunkier if you want the procedure factory to contain literal sin(x*y) instead of the reference name expr.

A variant on G2 above:

@sand15 That procedure F uses unapply to get an operator form for the integrand.

The result of
  unapply(expression_in_x, x)
is an operator (procedure).  As such, it no longer has a direct dependency on the name x in particular, or on any other name.

And so those commands like Int, fsolve, plot, etc, would take a mere range as second argument, as opposed to how they'd take name=range to accompany an expression.

In the case of Int, the numeric integration is done with respect to the arguments that will be passed to the operator-form integrand. In the case of fsolve, the root-finding is done with respect to that. Etc.

Forgive me if the following is too long-winded. (It's the kind of thing I'd put early in a book on Maple, if I were to write one.) You might well already know all this.

Download expression_vs_operator_form.mw

@Hullzie16 The optionsimplicit=[gridrefine=1] gets passed by plots:-inequal to plots:-implicitplot (which the former uses, internally). You could experiment with other implicitplot options there, to get a better resolution in the inequal plot. (That can make it take longer, say.)

The various digits, epsilon, and method options are used in the Int calls. Those get passed to the numeric integrator when those Int calls get evalf'd. Here it's tricky because you have nested Int calls.

In such situations the inner numeric integration often needs a more accurate result (smaller epsilon) to that the outer numeric integration can get its convergence for each evaluation of its own, outer integration.

The plotting command (or implicitplot) may be able to handle a coarser epsilon for the outer integration, which may make that faster. That may also be necessary to get the outer integration to converge (for the various values of the plotting variables).

Balancing all that, and the methods, can be tricky. Unfortunate choices may make it take a very long time, or not work at all.

@NIMA112 This is not super fast, and some computational effort is duplicated.

But it allows you to find the ranges for h[1] and h[2], and then merge them into a common range from which shared contour values may be constructed.

h1_h2_ok_ac.mw

@C_R My original choice of op was unnecessarily confusing, sorry.

I've edited my Answer. Hopefully it's more clear now.

Let me know if you still have queries on it that Preben has not explained.

Best of the season.

ps. Other people (eg. Joe Riel) have pointed out in the past that there are situations in which it would be programmatically convenient if there were a named,  stock "identity" function, ie. a stock name for x->x .

If it were called, say, ident then all I'd have to do is set `print/ident` to Typesetting:-Typeset and then utilize it here will calls like %ident(1) .

@Scot Gould It might also be possible to construct an appliable module which accepts two arguments. Through Explore, its first argument could be hooked up to a Slider while its second argument could be hooked up to a listbox/combobox/radio-button.

The Slider might provide 1 or 2 digits (to the right of the decimal mark) at a time, while the other control might specify the current scale. The appliable module could store the running total, eg. a+b/100+c/10000 or what have you.

That might be done reasonably cleanly, and tightly, hopefully easy to use.

I suppose that a neat and tidy set of such dymanic controls might (ideally? one day) customize any Explore Slider, with a hidden State component storing the running values of each scale.

@sand15 Yes, this is the very kind of thing that I meant when I referred to splitting a parameter.

I suspect that if the Sliders were long enough (and maybe snap-to-ticks) then one might be able to usefully & reliably get two powers of 10 (ie. two digits) from each.

So that'd allow for pinpoint specification of 10e-4 precision with just a pair, and might be quite nice.

I have another idea, but it'll have to wait until tonight...

In your old Maple 13 you might also look at something like, say,

   numtheory[sigma](n) - n