Use eliminate for a symbolic solution. For your system above, eliminate can solve for any three of the variables in terms of the fourth.

for V in combinat:-choose({x1,x2,y1,y2}, 3) do eliminate({eq||(1..3)}, V) end do;

You are mixing concatenation with indexing. If you want to consistently use indexing, then do

f_[0]:= x-> r(x):
for i to 10 do
     f_[i]:= subs(_I= i, x-> r(x)*f_[_I-1](r(x)))
end do:

(There is no need for 1 to be a special case.)

If you want to consistently use concatenation, then do

f_0:= x-> r(x):
for i to 10 do
     f_||i:= subs(_I= i, x-> r(x)*f_||(_I-1)(r(x))) 
end do:

Note that a name created by the concatenation operator || cannot be a local.

With a parameterized procedure definition (that is, a procedure created with a variable (the parameter) for which you wish to supply a value at the time the procedure is created), you need to use subs to supply the parameter's value. If you don't, the parameter remains symbolic even if it had a value at the time the procedure was created (unless the parameter is a module name (any other exceptions? Acer?)).

An alternative to subs is to use unapply to turn fully evaluated expressions into procedures. It's not really appropriate in this context (if r is complicated the procedures will be horrendously long), but it's done like this:

f_[0]:= x-> r(x):
for i to 10 do
     f_[i]:= unapply(r(x)*f_[i-1](r(x)), x)
end do:

Your procedures f_[i] (in the indexed form) can be defined as a single, simple procedure without the need for a loop or recursion. In the following, pay attention to how the index with which the procedure is called can be accessed within the procedure by using op(procname).

f_:= proc(x)
local i:= op(procname), k;
     if not procname::indexed or not i::nonnegint then
          error "This procedure needs a nonnegative integer index, e.g., f_[3](x)."
     end if;
     mul((r@@k)(x), k= 1..i+1)
end proc:

At this point, you might as well make i a declared procedure parameter rather than an index. Then the whole thing is even simpler---a one-liner---because no explicit error checking is needed:

f_:= (x, i::nonnegint)-> mul((r@@_k)(x), _k= 1..i+1):

This last form is what I would use.

Your error is in this passage:

Comet := array([1, 2]); if j = 10 then Colors := [aquamarine]; Linestyle := [longdash]; Comet[1] := spacecurve([subs(E[j] = E, x[j]), subs(E[j] = E, y[j]), subs(E[j] = E, z[j])], E = 0 .. 2*Pi, color = Colors[j], linestyle = Linestyle[j]) end if;

You redefine Color to be a list with one element, and then you ask for Colors[j] where j=10. The lists Colors and Linestyles already have 13 elements at this point in the code. There's no need for you to redefine them. The equivalent error occurs in the Oort passage immediately below the above passage. Simply remove the redefinitions of Color and Linestyle from these passages.

There's still an error after that; hopefully you will find it. If not, ask.

It would be helpful if you uploaded your work as a worksheet. If you do use plaintext, at least remove the > from the beginnings of the lines. If you don't, I have to remove every one of those by hand before I can run your code.

The LinearAlgebra:-SubMatrix command is superfluous. The same results can be obtained with indexing alone:

A:= <1,2,3,4; 5,6,7,8; 9,0,1,2>;

A[[1,2],[2,4]];

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.