An expression such a 'InverseErf/erf' needs to have both name quotes and unevaluation quotes, with the unevaluation quotes on the outside: '`InverseErf/erf`'(x). Without the name quotes, the expression is interpreted as an unevaluated division. This is not necessarily your only error; it is just the error that I saw immediately.

When you are inside a procedure F, and you want to return 'F'(x), it is preferable to return 'procname'(x).

It is not the same error. You are much closer to the answer now.

The current problem is that points1[2..6] contain a huge number of integrals that cannot be numerically evaluated because M has no value. In the commands eval(mH, {r= R, epsilon= eps1[i]})  you also need to specify a value for M.

Above the line b:= 2*Pi, add a line a:= whatever. Change the 2*Pi to whatever you want. Change every other occurence of b to (b-a).

applyrule:

Your second applyrule does not work because it gets caught in an infinite loop. It is very easy for applyrule to get into an infinite loop: All that you need is for the right side of the rule to match the pattern. Here's a quote from ?applyrule :

It applies the rules until no rule can be applied any more.

That being said, I don't know why your first applyrule does not also get into an infinite loop. But I can modify it so that it works:

collectExp:=proc(v_exp, v_collect)
local c1, c2:
     value(applyrule(
          c1::algebraic*exp(c2::algebraic)= c1*%exp(%collect(c2, v_collect)),
          v_exp
     ));
end;

The %s change exp and collect to unevaluated forms, and the value evaluates them (and changes %exp back to exp). I need to % on the collect so that it actually does the collect (don't know why), and I need the % on exp to prevent an infinite loop (because %exp will not match the pattern).

evalindets:

I usually find applyrule too finicky for my purposes, and I use evalindets (or subsindets) instead. If you want to apply the rule to arguments of exp only, do this:

evalindets(a1, exp(algebraic), ex-> map(collect, ex, [x,t]));

That's identical to

evalindets(a1, exp(algebraic), curry(map,collect), [x,t]);

If you want to apply the rule to all function arguments, just change exp(algebraic) to function.

Assigning values to variables upon which assumptions have been made doesn't work. Remove all the assignments, M:= 1, ..., omega:= 1. Then do a restart. Then reexecute up to the point where you have the integral. Then enter the integral using eval:

mean_integral:= 2*Pi*(Int(
      eval(simplify(mean_square)*area_element,
           [M = 1, r = 3, epsilon = 1, alpha = 1, beta = 1, omega = 1]
      ),
      theta = 0 .. Pi, method = _Gquad
));
evalf(%);
                        65.6364799988902

First, do a restart. Then add the coordinate names [x,y] as the third argument to both the circle and line commands to avoid the annoying dialogs which ask you to name the axes.

If I enter x and y in the dialogs, then your code above works for me (in Maple 17.02), and I get these points:

coordinates(F1);
              [29.9701216287282, 223.283461472248]
coordinates(F2);
              [31.9486961310283, 223.278449723452]

fsolve(eval(temp, x= 2.2), Omega= 0..infinity);
                       0.857190359134816

You can even make a function out of it and plot it:
F:= xi-> fsolve(eval(temp, x= xi), Omega= 0..infinity):
plot(F, 0..4);

I can't find online documentation for normalcdf, so I'll try to guess what it does. Does this do the job?

normalcdf:= proc(x1,x2,mu,sigma)
local t, dist:= Statistics:-CDF(Normal(mu, sigma), t);
     eval(dist, t= x2) - eval(dist, t= x1)
end proc;

You'll have to show me an example of arrow not working with CrossProduct, because it works for me.

Download arrow.mw

Addressing your other question, about which arrow will be used: It will be the one most recently loaded with with. If there's any ambiguity, you should use the module prefix, either plots:-arrow or plottools:-arrow. You should do this in your procedure; you should not rely on with within a procedure. Or, instead, you may use a ?uses clause within a procedure.

Here is a not-fully-polished implementation of Horner's method. It will find the non-multiple real roots of a polynomial with real coefficients. It will likely fail if it encounters a multiple root. Its output is the roots followed by the remaining unfactored polynomial.

Suggestions for improvement are most welcome.

Horner:= proc(P::polynom)
local
    p:= expand(P),
    x:= indets(p, And(name, Not(constant))),
    oldDigits:= Digits,
    roots:= table(),
    R, d
;    
    # Check that p is univariate.
     if nops(x) > 1 then
          error "Multiple names found: ", x
     else
          x:= x[]
     end if;
     
     Digits:= Digits+2+ilog10(Digits);

     # Get upper bound of abs val of roots by Rouche's theorem
     p:= evalf(collect(p,x));
     p:= expand(p/lcoeff(p));
     R:= 1+max(abs~([coeffs(p)]));
    
     for d to degree(p) do
         try
            R:= Student:-NumericalAnalysis:-Newton(
                 p, x= R, tolerance= 10^(-oldDigits-1), maxiterations= 2*Digits
            )
        catch:
            userinfo(1, Horner, StringTools:-FormatMessage(lastexception[2..-1]));
            break
        end try;
        p:= quo(p, x-R, x); # Reduce degree of p using the root
        roots[d]:= R        
    end do;
    
    evalf[oldDigits](sort(convert(roots, list)))[], evalf[oldDigits](p)
end proc:  

Example:

p:= `*`('randpoly(x, degree= 3)' $ 3);

infolevel[Horner]:= 1:
Horner(p);
Horner: maximum number of iterations (36) exceeded

fsolve(p);
 

Suppose that your Array is named A, and that its columns are indexed starting with 1 (the default). Then to plot the first vs. second columns

plot(A[.., [1,2]]); or plot(A[.., [2,1]]);

and to plot the first vs. third columns

plot(A[.., [1,3]]); or plot(A[.., [3,1]]);

exports(foo, static);

See ?exports .

There should be separate pages for each object and each object-method combination. Redundacy is acceptable in the "Description" sections, but try to vary the examples. Do not put things on the same page just because they are the same name overloaded.

This should be a Post, not a Question, because the answers are open-ended rather than definitive. Also, Posts tend to get more responses than Questions.

Yes, your commands are correct. But the answers returned by implicitdiff are more complicated than they need to be. After implicitdiff, use factor(simplify(%)).

See the examples at ?maximize or ?minimize .