You are ignoring the warning messages that Maple is giving you about variables being implicitly declared local. You didn't even mention them in your post, as if you thought those warnings were innocuous and completely irrelevant. Your arrays need to be declared as globals in procedure code. So, it should look like

code:= proc()
global U,R,m, ... and all the other Arrays ;

You can specify a domain restriction to fsolve. Change your fsolve command from this (your original):

poleR:= (M,g)-> fsolve(unapply(Denom(s, M, g), s));

to this:

poleR:= proc(M, g, R::range:= NULL)
     fsolve(unapply(Denom(s, M, g), s), complex, R)
end proc:

Then you can optionally pass a domain restriction for the solution:

poleR(3, .2); #No domain restriction.

If you're looking for the unique real root near 9 that you know exists, do

poleR(3, .2, 8..10);

     8.998575612

I used Maple 16. I can't test that in Maple 13 for you. Let me know if it works for you. I don't know if Maple 13 had optional positional parameters, which is what R is in poleR. If it doesn't work, I can easily work out something else for you.

Use this Monte Carlo algorithm: Pick four points on the curve at random. Use the result of passing those four points to geom3d:-AreCoplanar.

There are several steps needed to make this work.

1. You need the output from RandomTools:-Generate to be a procedure that generates a random number rather than being just a single random number. You do this by including the option makeproc.

Ra:= RandomTools:-Generate(distribution(Uniform(0.001, 0.02)), makeproc);

2. You need c to be a procedure which returns a random number when it is passed a numeric argument and which returns unevaluated when passed a symbolic argument. You do that like this:

c:= t-> `if`(t::numeric, Ra(), 'procname'(args));

Note the two different types of quotes in that expression. Don't mix them up and don't omit them.

3. You need to refer to c as c(t) in the differential equation:

eq:= diff(X(t),t)= 1-f-c(t)*X(t);

4. You need to tell dsolve that c is a "known function" by including the option known= [c].

5. You need to defeat the dynamic step sizing used by modern IVP solvers because the randomness will cause it to use very small sizes and it will thus take a very long time to finish, if ever. You do this by choosing a classical IVP solver. I chose the most basic one, foreuler (forward Euler), but others should work. You can't use the range option with the classical solvers.

sol:= dsolve({eq, init}, {X(t)}, known= [c], numeric, method= classical[foreuler]);

Putting that all together, we have

Ra:= RandomTools:-Generate(distribution(Uniform(0.001, 0.02)), makeproc);
c:= t-> `if`(t::numeric, Ra(), 'procname'(args));
f:= .1;
eq:= diff(X(t),t)= 1-f-c(t)*X(t);
init := X(0) = 100;
sol:= dsolve({eq, init}, {X(t)}, known= [c], numeric, method= classical[foreuler]);
plots:-odeplot(sol, [[t, X(t)]], t = 0 .. 100);

There are at least three ways to do this.

type(expr, freeof(x));

     false

depends(expr, x);

     true

type(expr, dependent(x));

     true

Note, however, that if the expression contains x, but only as a bound variable, then freeof, depends, and dependent, will return true, false, false, respectively. If you do not care about the distinction between free and bound variables (i.e., you just want to know whether expr has any x at all), then use

has(expr, x);

(For those readers who may not have studied formal logic, an example of a bound variable is the x in any of the following:

Int(f(x), x= a..b);
Sum(f(x), x= a..b);
Eval(f(x), x= a);

etc.  It is a variable that would not appear in the evaluated result of the expression were it possible to fully evaluate it. Another way to think of it is that the variable is bound to the expression if it would not make any sense to assign it a value outside of the expression.)

I think that Matlab gives an answer because it is not being as cautious as Maple is trying to be. The trick to getting Maple to do this integral numerically is to get Maple to "throw caution to the wind." In this case, you do that by making the input to the Ints numeric procedures rather than expressions. Then it can't "look inside" the procedures, so it doesn't "think too hard" about them, but rather proceeds with its most basic numeric integration algorithms.

Another trick that helps here is to reduce the requested precision by using the epsilon option to Int.

[Edit: I originally had 11+1/2 where the code below has 5+1/2. That was a mistake.]

restart:
Digits:= 15:
f1:= (1-exp(-(5+1/2)/cos(y)))*sin(y):
yrange:= 0..arctan(300*cos(x)+sqrt((12+1/4)-90000*sin(x)^2)):
xrange:= 0..Pi:
f1y:= unapply(f1, y):
yrangex:= unapply(yrange, x):
evalf(Int(x-> evalf(Int(f1y, yrangex(x), epsilon= .5e-12)), xrange, epsilon= .5e-9));

The true value of the imaginary part is presumably 0, and thus it can be ignored as roundoff error. The real part should only be trusted to nine digits. You can increase the precision by increasing the value of Digits and decreasing the values of the epsilons.

Another way to get Maple to "throw caution to the wind" is to use option continuous on a symbolic integration of the inner integral. To forces Maple to apply the fundamental theorem of calculus even though it may not be able to verify continuity. I believe this is akin to Axel's way.

evalf(Int(unapply(int(f1, y= yrange, continuous), x), xrange, epsilon= .5e-12));

Note that the real part is the same for 13 digits, which is the precision specified by epsilon.

Here's how to do it. I plotted it for the third value of c0 and the third value of m. I hope that you see how you can do it for any combination of c0 and m.

X:= (c0,m)->
     c->
          sqrt(2/(m*c0^(m-2))) * (c/c0)^(1-m) * sqrt((c/c0)^m-1) *
          hypergeom([1,1-1/m], [3/2], 1-(c/c0)^(-m)) / sqrt(2)
:
C0:= [0,0.277164,0.340257,0.408617,0.581345,0.649268]:
M:= [0,1/2,1/3,3/4,2,3]:
plot([X(C0[3],M[3])(c), c, c= 0..1], view= [0..1, 0..1]);

The writedata will work if you use default as the fileID. But it is easier and more flexible to just use printf.

printf("%6.2e", 0.3*10^(-5));

     3.00e-06

If you insist that you want the leading digit to be 0, it can be done, but it is significantly more work. Here's a procedure:

MyPrint:= proc(x::float)
local mm,m,e;
     (m,e):= op(x);
     while irem(m,10,'mm')=0 do m:= mm; e:= e+1 end do;
     printf("%s0.%de%d", `if`(x<0, "-", ""), CopySign(m,1), e+length(m))
end proc;

MyPrint(.3*10^(-5));

     0.3e-5

You seem reluctant to simply enter the integral

int(c*g(c), c= a..b);

Why are you reluctant?

Yes, you can use Statistics:-ExpectedValue for an arbitrary distribution. But what's the point of doing it that way?

Suppose my distribution is g(x) = x/2 defined on [0,2]. I can do

Omega:= Statistics:-Distribution(PDF= (x-> x/2), Support= 0..2):
Statistics:-ExpectedValue(Omega);

or I can simply do

int(x*x/2, x= 0..2);

and get the same answer.

Is your problem that you have no closed form for the PDF? Then the Statistics package may be able to help because there are many alternative ways to specify a distribution. What information do you have about the distribution?

Your odetest result is {0, expression1, expression2}. Your input was four differential equations. In Maple a sequence of expressions surrounded by curly braces is a set. Sets don't have repeated elements. Two of your differential equations matched perfectly to their equations, producing 0. Those two zeroes are replaced by one zero because it's a set. To avoid this confusion, use lists instead of sets. Practically, that means replacing all your curly braces { } with square brackets [ ].

The expression1 and expression2 may well reduce to 0 also, but Maple was not able to prove it. Substitute numeric values for as many parameters as you can and see if you can reduce these to 0. I wanted you to post your FULL CODE so that I could try to reduce these to 0, or prove that they aren't 0. So, if you want more help, please post that.

Presumably, you intended for your procedure to be like this:

Naief:= proc(A::integer, B::integer, p::posint)
local x, a:= A mod p, b:= B mod p;
     for x to p-1 do
          if a^x mod p = b then return x end if
     end do;
     print("No solution.")
end proc:

Naief(3,7,11);

     "No solution."

Naief(7,3,11);

     4

Here is some Maple coding advice for you:

There is no disp; what made you think that there is? It's either print or return. In the vast majority of cases, you'll want return instead of print, which is mostly for auxiliary information. See also userinfo. Of course, return forces an exit from the procedure, and print does not.

If you want a null operation, the syntax allows you to leave the space blank; there's no need for something like x=x. In other words, the syntax allows then elif with nothing between them. If you want to immediately go to the next value of x, use next.

Don't check both a condition and its opposite; it is wasteful. A thing either is or isn't equal to 0. So, you should use else instead of elif.

Did you have any trouble reverting the substitution? A simple

eval(expr, Sw= Sum(w[k], k= 1..n);

should do it.

If the sum cannot be simplified or be put into closed form, then you should use Sum instead of sum. If you use sum, then Maple wastes time trying to put the sum into closed form everytime the expression is evaluated.

This is an Answer to the amended question in your Reply to your original Question, which is a much more complicated problem than your original Question.

evalindets(B, And(`*`, satisfies(P-> {1/T[1],1/V[2]} subset {op(P)})), P-> P*T[1]*V[2]);

(It looks a lot better in a worksheet than posted here.)

Use the Select Execution Group keyboard shortcut Ctrl-Alt-Shift-E. It's true that that will only select one execution group at a time, but, c'mon, how many execution groups are we talking about moving here? If you have a huge number of execution groups, maybe you should merge them using Join Execution Groups F4.

Note how Maple represents the derivative of tan(x):

diff(tan(x), x);

Based on this, I guessed that your integral might work better if I replaced sec(x)^2 with the above.

int(tan(x)^(n-2)*(1+tan(x)^2), x);

The above answer is returned nearly instantaneously.

I was quite surprised that the following not-quite-simplified answer was also returned nearly instantaneously. I had expected that the integration would be distributed over the sum, leading to an unevaluated answer for powers of tangent.

int(tan(x)^(n-2)+tan(x)^n, x);

Based on the result of my first integral above, I guessed that if I redefined the derivative of tangent, then your integral would evaluate easily. That guess was correct.

eval(`diff/tan`);

restart: #Need to clear diff's remember table.
unprotect(`diff/tan`);
`diff/tan`:= proc(z,t,$) `tools/df_dt`('tan(z)', z-> sec(z)^2, t) end proc:
diff(tan(x), x);

int(tan(x)^(n-2)*sec(x)^2, x);

So this explains why the derivative-divides rule (see Alejandro's Answer) failed to be applied to your original integral.