Your problem 1 can be solved by a single Maple command:

Interp([0,1,6], [1,5,2], x) mod 7;

It's entered just like your problem from earlier today.

PDE:= diff(r*diff(T(r,t),r),r)/r = diff(T(r,t),t)/beta:
pdsolve(%, build);

Here's a solution that makes use of parallel processors.

PowerSet:= proc(S::set)
local P,x,s;
     P:= {};
     for x in S do  P:= P, Threads:-Seq({s[], x}, s= {P})  end do;
     {P}
end proc;

pdsolve(diff(y(x,t),x$2)+diff(y(x,t),x)/x=diff(y(x,t),t)/alpha);

Here is why your subs did not work. Consider the expression x+2+sin(x+2) as a tree (in the mathematical sense of tree). The root is the whole expression. The first level of branching divides into three nodes: x, 2, and sin(x+2). The x and the 2 are terminal (leaf) nodes. The sin(x+2) splits into two nodes: the sin and the x+2. The sin is a leaf, and the x+2 splits into two leaf nodes. (Pause now to sketch the tree.) The expression which gets substituted for, i.e., the left side of the = in the subs command, has to be a node in that tree. The x+2 only appears as a single node when it is the argument of sin. The other x and 2 get split immediately from the root.

Does that makes sense?

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);