The following paragraph is from the help page ?dsolve,details:

Integrals appearing in answers returned by dsolve are expressed using the inert Int and Intat (not int or intat). These integrals appear when int is not able to calculate them or when it appears to be convenient not to evaluate them during the solving process. You can request the evaluation of these integrals using the value command.

So, try applying the value command to dsolve's results. If you still get an integral, then Maple simply can't do the integral. If you have parameters, sometimes making assumptions on the parameters will make the integral doable. See ?assuming.

If you have further questions about this, please post your differential equation and any initial or boundary conditions.

The $ is the repetition operator (it also has a mode where it acts like seq). It is a basic operator of Maple, not something special to permute. For example, [1$4] returns [1, 1, 1, 1]. In many ways, $ is to seq as sum is to add and as product is to mul. For example, $ will take a symbolic second argument. (It's a shame that seq and $ are not switched: The former is more commonly used.)

To get all 2^n sequences of 0s and 1s of length n, use

combinat:-permute([1$n, 0$n], n);

While the $ does accept symbolic arguments, permute does not. To run your command, you'll need to supply a nonnegative integer for n.

Pi is not considered numeric. So, you have to do something like this

for x[2] from 15/180 to 75/180 by 5/180 do    
     a[i,j]:= Pi*x[2];
     j:= j+1
end do;

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