When you need scan through an entire rtable (any Vector, Matrix, or Array), probably the best command to use is rtable_scanblock. This command is builtin and extremely fast. The following worksheet contains code to solve your reverse-index problem and a comprehensive collection of convenient access methods. It addresses all the issues discussed by the other respondents. I doubt that there can be anything faster.

This code builds the reverse index once, and quickly, so that all subsequent inquiries are essentially instantaneous. It also handles the lookup of individual entries, subsets of entries, non-existent entries, or any combination of those.

RtableReverseIndex.mw

MaplePrimes is refusing to display the worksheet inline, so you'll need to download it to get all the details. The main block of code in the worksheet uses a feature new to Maple 2018. If you're using an earlier version, I can make a minor change (it'd only take 2 minutes) so that it'll work in earlier versions. Let me know.

Here's the worksheet in plaintext with the output removed:

restart:
This procedure contains everything needed to make a reusable reverse index for any rtable (any Vector, Matrix, or Array) of any number of dimensions.
RtableRI:= proc(A::rtable, $)
local T:= table(), v;
   rtable_scanblock(A, [], (v,j)-> (T[v][j]:= ())); #Build reverse index.
   for v in indices(T, 'nolist') do #for each entry...
      #...count indices; make the entry's indices a list; store both in a vector:  
      T[v]:= (J-> <nops(J)|J>)([indices(T[v], 'nolist')])     
   od;
   eval(T)
end proc:

The above procedure uses a coding feature (embedded assignment) that was only made available in Maple 2018. If you're using an earlier Maple version, use the following:
RtableRI:= proc(A::rtable, $)
local T:= table(), v;
   rtable_scanblock(A, [], proc(v,j) T[v][j]:= () end proc);earlier 
   for v in indices(T, 'nolist') do
      T[v]:= (J-> <nops(J)|J>)([indices(T[v], 'nolist')])     
   od;
   eval(T)
end proc:

It works like this....
A:= <a,a; b,c>;
T:= RtableRI(A);
Lookup entry a:
T[a];
It occurs 2 times and at the given indices. 

A tabular presentation of the entire table can be made thus:
M:= (T-> <<lhs~(T)> | <rhs~(T)[]>>)([entries(T, 'pairs')]);
And a DataFrame can used for a slightly more-elegant presentation:
DF:= DataFrame(op~(M[..,2..]), rows= M[..,1], columns= [Count, Indices]);
All of the above, plus some additional convenient access methods to the data, are collected together in the following module:
RtableReverseIndex:= module()
description
   "Builds a reverse index of any rtable so that subsequent lookup of entries "
   "is essentially instantaneous. Lookup results are returned as DataFrames."
;
option 
   `Author: Carl Love <carl.j.love@gmail.com> 23-Aug-2018`, 
   object
;
local 
   #just a convenient abbreviation:
   Inds::static:= proc(T::table) option inline; indices(T, 'nolist') end proc,

   #common output row for non-existent entries, if they're searched for:
   ZRow::static:= rtable(1..2, {1= 0, 2= ()}),

   ###### some stuff for formatting DataFrames ##########################################
   Cols::static:= 'columns'= ['Count', 'Indices'],                                      #
   Names::static:= proc(J::{list,set}) option inline; convert~([J[]], 'name') end proc, #
   Rows::static:= (J::{list,set})-> 'rows'= Names([J[]]),                               #
   ModulePrint::static:= (A)-> A:-DF, #Prettyprint raw objects as their DataFrame.      #
   ###### end of formatting stuff #######################################################

   #main access point: the new-object-constructing procedure:
   ModuleApply::static:= proc(A::rtable, $)::RtableReverseIndex; 
   local New:= Object(RtableReverseIndex);
      New:-DF:= DataFrame(
         <entries((New:-Tab:= RtableRI(A)), 'nolist', 'indexorder')>, 
         Rows((New:-Ent:= {Inds(New:-Tab)})), Cols
      ); 
      New
   end proc
;
export
   ###### the dynamic part of the object #
   Tab::table,    #the master index      #
   DF::DataFrame, #the master DataFrame  #
   Ent::set,      #the rtable's entries  #
   ###### end of dynamic part ############

   #the index-creating procedure:
   RtableRI::static:= proc(A::rtable, $)::table; 
   local T:= table(), v;
      rtable_scanblock(A, [], (v,j)-> (T[v][j]:= ())); #Build index.
      #Count each entry's indices and convert them to a sequence.
      for v in Inds(T) do T[v]:= (()-> rtable(1..2, {1= nargs, 2= args}))(Inds(T[v])) od; 
      eval(T)
   end proc,

   #the index-accessing procedure: 
   `?[]`::static:= proc(A::RtableReverseIndex, J::list, $)
   local 
      J0:= {J[]} minus A:-Ent,
      R:= `if`(
         J0 = {}, #all requested entries exist
         A:-DF,
         #Append a row for each non-existent entry:  
         Append(A:-DF, DataFrame(`<,>`(ZRow $ nops(J0)), Rows(J0), Cols))
      )[Names(`if`(J=[], A:-Ent, J)), ..]
  ;
      `if`(nops(J)=1, [R[1,1], [R[1,2]]][], R) 
   end proc 
;
end module: 
      
Some examples:
interface(rtablesize= 30):
randomize(1): #just to stabilize these examples
A:= LinearAlgebra:-RandomMatrix((5,5), generator= rand(-5..5));
ARI:= RtableReverseIndex(A);
You can request the sub-table for any subset of entries, in any order, including entries that don't exist:
ARI[4,-1,7,1];
If you request a single entry, then just the count and the list of indices is returned, without the surrounding DataFrame:
ARI[4];
That applies to non-existent entries also:
ARI[7];
Test on a large example:
A:= LinearAlgebra:-RandomMatrix((2^10)$2, generator= rand(-5..5));
So, over a million entries. 
ARI:= CodeTools:-Usage(RtableReverseIndex(A)):
Once the index is created, all subsequent lookups are essentially instantaneous. Here, I extract the column of counts:
CodeTools:-Usage(
   ARI[][Count]
);
And here I extract (without displaying) the list of all indices for entry 0:
ZeroIndices:= CodeTools:-Usage(
   ARI[0][2]
):
nops(ZeroIndices);
A small sample of those:
ZeroIndices[100..150];

You can simply integrate your expression. There's no need to simplify anything or assume that anything is real. The following returns 2*Pi:

int(conjugate(exp(I*phi))*exp(I*phi), phi= 0..2*Pi);

Also, the product of a complex function and its conjugate is the square of its absolute value (regardless of whether the parameters are real). So, the following also returns 2*Pi:

int(abs(exp(I*phi))^2, phi= 0..2*Pi);

You are trying to compute the expected value (E(.)) of a linear combination of r.v.s. As is very well known, expected value is a linear operator. So the answer is 0.3*E(p1) + 0.7*E(p2). Now, you can get this with Statistics:-Mean, as has been pointed out, or you can even use a pocket calculator: 0.3*1/(1+100) + 0.7*1/(1+50).

Maple's Statistics package is not very good at finding the PDFs of combinations of r.v.s, even linear combinations. I suspect that this has something to do with difficulty simplifying combinations of the piecewise expressions that are used for the PDFs of distributions whose support is not the whole real line.

For array input: As far as I can tell, it's not directly possible. But you can put your array data through CurveFitting:-Spline to turn it into a suitably smooth piecewise polynomial that interpolates the data.

For array output: The easiest thing is to make a plot and extract the data matrix from it. For most plots P, the data matrix is op([1,1], P).

If I define w(x) = int(y*u(y)^3, y= 0..x) and w(1) = J, then your IVP becomes the BVP system

eq1:= diff(w(x),x) = x*u(x)^3, w(0)=0, w(1)=J;
eq2:= diff(u(x),x) = exp(x) + (2*exp(3)-1)/9 + J, u(0)=1; 

Hopefully the derivation of that is obvious; if not, let me know. The system has a parameter, J. Maple will numerically solve BVPs with parameters provided that there's an additional BC for each parameter. As you can see, we have 3 BCs and the differential order is 2; so, we're good to go.

dsn:= dsolve({eq1, eq2}, numeric);
Error, (in dsolve/numeric/bvp) initial Newton iteration is not converging

Uh oh, we hit a roadblock. This is a very common error message from Maple's BVP solver. When you get this error (note the word initial), the first thing to try is a kind of bootstrapping called a continuation parameter. This is an additional parameter C that varies continuously from 0 to 1 that you place anywhere in the system (usually as a coefficient) such that when C=1, you get the original BVP, and when C=0, you have a relatively simpler BVP (see ?dsolve,numeric,BVP). It may take some guessing and a few trials to find a place to put C that works. (You can also put C in multiple places at once.)

In this case, what works is changing the BC w(1)=J to w(1)=J*C:

eq1:= diff(w(x),x) = x*u(x)^3, w(0)=0, w(1)=J*C;
eq2:= diff(u(x),x) = exp(x) + (2*exp(3)-1)/9 + J, u(0)=1;
dsn:= dsolve({eq1, eq2}, numeric, continuation= C);
plots:-odeplot(dsn, [x,u(x)], x= 0..1); 

It works, and the resulting plot and a deeper analysis of the numeric values shows complete agreement with Mariusz's symbolic solution.

If you need to continue the solution for u(x) beyond the interval 0..1, that can be done with a little extra work.

Okay, here is a procedure to do it:

VuReduce:= proc(e)
uses P= Physics;
local old:= P:-Expand(e), r;
   do
      r:= (P:-Expand@evalindets)(
         old,
         specfunc(P:-`*`),
         proc(`*`)
         local p,p1,U,V;
            if
               membertype({v[nonnegint], 'P:-`^`'(v[nonnegint], posint)}, `*`, 'p')
               and nops(`*`) > p
               and (U:= op((p1:= p+1), `*`))
                  ::{u[nonnegint], 'P:-`^`'(u[nonnegint], posint)}
            then
               V:= op(p,`*`);
               subsop(
                  p= Vu(
                     `if`(V::v[nonnegint], op(V), op([1,1],V)), # "a"
                     `if`(U::u[nonnegint], op(U), op([1,1],U)), # "b"
                     `if`(V::v[nonnegint], 1, op(2,V))          # "c"
                  ),
                  p1= `if`(U::u[nonnegint], NULL, P:-`^`(op(1,U), op(2,U)-1)),
                  `*`
               )
            else
               `*`
            fi
         end proc
      );
      if r = old then return r fi;
      old:= r
   od
end proc:

Call it via

VuReduce(e)

where e is any expression (or container of expressions) that you want to apply the transform to. Please test this. I made it so that it will  keep applying the transform as long as the expression keeps changing under it. Hopefully this won't lead to an infinite loop, but I can't know that without knowing the details of Vu. If it's a problem, I can easily include code to detect cycling. Note that I apply Physics:-Expand to both the original and transformed expressions.

If you have any trouble with this, please post a worksheet showing the erroneous output or error message and an explanation of what you think the result should've been.

I don't know if you're in the habit of using either 1D or 2D Input. I only use 1D. I haven't tested the above procedure entered as 2D Input. If this makes any difference, it would only affect the entry of the procedure itself. The input form of the expressions or the procedure's call won't matter. I'm not expecting there to be any differences, but there are many unfortunate surprises with how 2D parses the code.

I've never worked with the Physics package before, so there may be some easier way to do this that I missed. Nonetheless, the above should work with reasonable efficiency.

You could do it like this:

Trunc:= proc(F::{procedure, `+`}, odr::posint:= 2, v::list(name):= [x,y,z])
local f:= `if`(F::procedure, F(v[]), F);
   if not f::`+` then
      error "Can't truncate 1-term expression"
   else
      select(q-> degree(q,v) <= odr, f)
   fi
end proc:

As hinted at by Tom, the only purpose of ArrayTools:-Alias is to create a reshaped, re-indexed, resized, and/or re-ordered secondary reference to an rtable (an Array, Vector, or Matrix) or to some subpart of it. If you apply Alias directly to a Matrix(...) command, then there's no primary reference to that Matrix, so there can be no point to creating a secondary reference (although it will allow you to do so). A primary reference would be created by first assigning the result of Matrix(...) to a variable.

You have not declared Wn to be an Array. When an assignment is made to an indexed undeclared variable, a table (but not an rtable) is created. But table indexing works differently than Array indexing in that index ranges such as 0..1 are not respected. In your case, Wn ends up with only 3 indices---[0..1, 0..1, 1], [0..1, 0..1, 2], and [0..1, 0..1, 3]---so the ranges are literally part of each individual index. You can inspect the end result of your loop with eval(Wn).

Representing  (1 - 1e-18) requires 18 significant digits, whereas representing 1e-18 requires only 1 significant digit. So I find the following way much easier than all the others presented in this thread, requiring no increase in Digits from its default. I apply the following obvious identity:

Quantile(Normal(mu,sigma), 1-x) = mu - sigma*Quantile(Normal(0,1), x).

Thus

Quantile(Normal(0,1), 1 - 1e-18) = -Statistics:-Quantile(Normal(0,1), 1e-18);
      8.75729034878321

The following procedure will apply the rule only to those lns that occur as a result of the integration:

REALINT:= proc(f::algebraic)
local OrigLn:= indets(f, specfunc(ln));
   evalindets(
      int(args), 
      specfunc(ln), 
      L-> `if`(L in OrigLn, L, ln(abs(op(L))))
   )
end proc:

Example:
REALINT(1/(x*ln(x)), x);

If you're willing to have Maple tacitly make the necessary assumptions that'll get you that most simplification, you can use simplify(..., symbolic). Although I've never seen this occur in a practical situation, it is possible to construct expressions where the necessary assumptions are contradictory, and thus the simplification is not valid for any possible valuations.  See ?simplify,symbolic.

Manipulating square roots is difficult. Here's one tip: Almost all sqrt(...) are converted to (...)^(1/2), so searching expressions for sqrt won't work.

Here's a procedure to integrate recursively until either there's a term with no derivative or we get stuck with an unevaluated integral:

MyOdeInt:= proc(ODE::algebraic, Z::name:= NULL)
local 
   ode:= frontend(expand, [ODE]), #Distribute, but don't expand functions.
   #Find the integration variable:
   z:= `if`(Z=NULL, map2(op, -1, indets(ode, 'specfunc'(diff)))[], Z),
   r 
;
    if ode::`+` then #There's more than 1 term.
       if remove(hasfun, ode, diff) <> 0 then return ode fi 
    #single-term case:
    elif not hasfun(ode, diff) then return ode 
    fi;
    r:= int(ode, z);
    `if`(hasfun(r, int), r, thisproc(r))
end proc:

And you'd use it:

MyOdeInt(odetemp);

I tested on higher orders:

MyOdeInt(diff(odetemp, z$9));

which returned the same thing as the first command.

At the beginning of your worksheet, do

with(Units); with(Standard);
UseUnit('J/(mol*K)'):
UseUnit('J/(kg*K)'):

Then computations whose results can be expressed in those units will be. Computations include arithmetic, evalf, and combine(..., units). The raw extraction of R from ScientificConstants will stay the same, but you can force it to the preferred units with simply combine(..., units) rather than needing convert.

HFloat stands for "hardware float". These are the 64-bit representations of floating-point (or decimal) numbers that are nearly universally used by computer chips. (They're called "double precision" in many languages.) Included are representations of infinity, -infinity, and undefined.

For most purposes, you can think of HFloat(infinity) as simply infinity. If it appears as a result of integration, it indicates some type of divergence issue. Perhaps the integral is truly divergent. If it appears multiplied by some coefficient expression C, and C could possibly be 0, then you may have an indeterminate form, which'll require deeper analysis. There's no simple way to deal with this. It requires case-by-case analysis.