@Axel Vogt 

I don't understand the purpose of your reformulation. Yes, it is interesting that Maple can perform that symbolic integration of the MeijerG expression in an instant. Perhaps that was your point. But the Int/diff expression in Question doesn't use the integral of the MeijerG expression; rather, it uses the integral of the square of the derivative of the MeijerG expression. This is much more complicated.

@Axel Vogt 

Here is my simplest "clean" interpretation of the Int/diff expression. By "clean", I mean that each of the three nested definite integrals has its own dummy variable of integration which is not used outside that integral:

H3p_Clean:= proc(ms::algebraic, b::numeric, eps::positive:= 0.5*10^(1-Digits))
local x:= indets(ms, And(name, Not(constant))), D, X, Y;
     if nops(x)>1 then
          error "1st argument must be an expression of at most 1 variable"
     end if;
     x:= x[];
     D:= diff(ms,x);
     evalf(
          subsindets(
               Int(
                    diff(
                         subs(x= Y, D)*Int((-b*X+1)*Int(D^2, x= 0..X), X= 1..Y),
                         Y
                    )*subs(x= Y, ms),
                    Y= 0..1
               ),
               specfunc(anything, Int),
               J-> Int(op(J), epsilon= eps)
          )
     )
end proc:

(I hope that my indentation, line break, and parenthesis style serves to elucidate rather than obfuscate the expression for you.)

In all cases that I've tried, this produces the same numeric results as the original "sloppy" version (posted in the Answer above). (By "sloppy", I mean that the same variable is used for all integrations and differentiations, and it's also used in the limits of the inner integrals.) But the weird thing is that in all cases the sloppy version ran significantly faster, usually by a factor of 2. So maybe the sloppy version is better.

I tried all of the following for the first parameter, ms: x, x^2, x^3, sin(x), cos(x), exp(x), exp(-x^2), sin(x)+cos(2*x). For each, I used both 0 and -1 for the second parameter, b. I also used the original MeijerG expression with b = .7. I varied the order in which the two integrals were done.

I have Digits set to 15, but I never used an epsilon smaller than .5e-10.

You could do an animation where each frame was one of a sequence of parallel planes intersecting the solid. In each plane, color would represent the value of the function. If you've ever seen an animation of a CAT scan or MRI, you'll know what I mean.

@frajola You wrote:

I recognize that it is a mathematical error to have the same variable in the limits of integration and on the integrand. Not only Maple does not care about it but also most researchers in the particular field that I am studying on.

It is very important that you check that Maple's interpretation of the integral is the same as that of the researchers in your field. Use a known simple case for ms; there's no need for any MeijerGs; you could use ms = x or x^2. Use b = 0 or -1. Compare Maple's answer with the known answer.

@Axel Vogt According to the fourth paragraph of Description at ?RootFinding,Isolate, the user may choose between the RegularChains Maple package or the RealSolving C library. You'd probably be interested in looking at the references at the end of ?RegularChains and ?RegularChains,RealTriangularize. The RegularChains package is massive, with several thousand procedures (of course, only a tiny percentage are documented).

@albert_heinle 

I don't know whether you'd consider this less hacky than specfunc, but you could do this

Ans:= fsolve(...);
if  Ans::function then
    #No real solutions found
else
      #Ans contains some real solution
end if;

It is intended that the above be done in a procedure. If done in top-level code, the use of eval(Ans,1) is recommended to avoid the re-evaluation of the fsolve.

Note that when fsolve (without option complex) returns unevaluated, it is not guaranteeing that there are no real solutions. It just means that it was unable to find any. On the other hand, RootFinding:-Isolate is attempting to give you a guaranteed answer. That's why it takes longer.

@acer 

Can you provide a simple example of an integration where assuming continuous makes a difference? It seems to differ from using continuous as an argument to int. Compare

int(floor(x), x= 0..X) assuming continuous;

with

int(floor(x), x= 0..X, continuous);

@Preben Alsholm 

The following plot of residuals, using fdiff to approximate the derivative, done in Maple 18, doesn't inpsire much confidence in dsolve's solution outside the interval 0.5..0.85. How does this turn out in Maple 2015?

eq:= diff(y(x),x) = cos(x) + y(y(x)-2):
res:= dsolve({eq, y(0)=1}, numeric, delaymax= 1, output= listprocedure, relerr= 10^(2-Digits)):
Y:= eval(y(x), res):
plot(fdiff(Y,[1],[x]) - subs(y= Y, rhs(eq)), x= 0..1);

I don't have much experience with fdiff, and I'm not absolutely sure that I'm using it correctly.

@NickB I edited the plaintext expressions into your original Question. The math editor in the MaplePrimes editor produces unreadable results. I don't know of any good use for it.

@maple fan 

You coded:

diff(f1(x,y),x)+diff(f1(x,y),y)+diff(f2(x,y),x)+diff(f2(x,y),y);
remove(has,%,indets(%,'diff'(anything,identical(x))));

The situation where the terms in their entirety are the derivatives is a very special case. In that case, using remove is equivalent to setting the derivatives to 0. But in the general case, using remove(has, ...) on an expression of type `+` would remove the entire term if the term had the derivative.

@maple fan 

If you want to make them equal to 0, use evalindets: Assuming that ex is your expression and that you want to set all derivatives with respect to x to 0, do

evalindets(ex, 'diff'(anything, identical(x)), 0);

Are you looking for a numeric or a symbolic solution to your differential equations? You have three nonlinear first-order ODEs each with a little over 100 terms. Maple will try, perhaps for a very very long time and using very much memory, to solve these symbolically. Or you can get a numeric solution in a second. Either way, I know how to correct the error.

I simplified the Boolean condition by using xor. Here it is:

`&<`:= proc(x::~Vector, y::~Vector)
local n:= op(1,x), i;
     for i to n while x[i] = y[i] do end do;
     i > n or not (i::even xor x[i] > y[i])
end proc;

Both this and the previous version are intended to be run in 1D input mode. The 2D input mode has different evaluation rules for Boolean expressions. This 1D-2D difference may have been fixed in Maple 2015; I can't check that right now.

If you need high-speed performance from this procedure, let me know. The speed can be improved at the expense of the flexibility by avoiding the Vector conversion and by setting the length to a constant.

@maple fan 

Well, what do you mean by "remove"? There's only two ways that I can think of to remove a subexpression from an algebraic expression. One is to set the subexpression to 0, and the other is to set the subexpression to 1. But there might be other ways that I'm not thinking of. Please show me an example of an expression before and after the derivatives have been removed.

@tomleslie 

In any Maple program, the sequence of random numbers will be the same unless you use the command randomize().

For 1000 trials, the probabilty of deviating from 667 by more than 32 is exactly

1 - add(binomial(1000,k)*(2/3)^k*(1/3)^(1000-k), k= 667-32..667+32);

or about 3%.

For two sets of 1000 trials, assuming properly generated random numbers, the probability of getting exactly the same number of successes is exactly

add((binomial(1000,k)*(2/3)^k*(1/3)^(1000-k))^2, k= 0..1000);

or about 2%.

The standard deviation of the underlying distribution is exactly

sqrt(1000*(2/3)*(1/3));

or about 15. Thus a deviation of 32 from 667 is slightly more than two standard deviations from the mean.