This is a bug in Maple's 2D Input. It misinterprets &xor as &x or. I suggest that you switch to 1D Input (aka Maple Input). Using xor instead or &xor will not (in general) give you correct results when using the Logic package! You need to use the & forms of all logical operators because the rest of Maple will consider them inert and only the Logic commands will process them.

2D Input is complete garbage with no redeeming value. However, if you insist on using it, you can work around this bug by using `&xor`(a,b) wherever you would've used a &xor b.

You have an eigenvalue with algebraic multiplicity 4. Depending on the numeric values of the parameters, the matrix may be defective---that is, have a eigenvalue for which the number of linearly independent eigenvectors is less than the algebraic multiplicity of the eigenvalue.

It may be possible to break down the cases into elaborate nested piecewise expressions. These'll specify the eigenvectors for all possible inequality ranges of the parameters. It may be possible to do the required polynomial^[1] algebra with solve(..., parametric...) (see ?solve,parametric), which is essentially the same thing as SolveTools:-Parametric (see ?SolveTools,Parametric). Then it may be possible to process each case with LinearAlgebra:-NullSpace to get the eigenvectors. If you decide that this is worth pursuing (and note where I've emphasized "may"), the first step is to specify allowed ranges for all parameters and to state whether they are real. This'll reduce the number of cases.

[Edit: The green part has been added; the red part deleted.]

It may also be possible to get a "generic" solution: one where it's implicitly assumed that all divide-by-zero^[2] conditions are met such that there are four linearly independent eigenvectors corresponding to the eigenvalue of multiplicity 4. This'll be easier than what I described in the previous paragraph, but we'll be assuming things that may not be true. That may not matter if you're going to substitute approximate decimal values for the parameters.

^[1]These commands are only implemented for polynomials. To apply them to rational functions requires that the denominators be cleared and have their zeroing conditions analyzed separately.

^[2]Row reduction of a matrix of polynomials requires division by polynomials. Careful case-wise analysis requires considering separately the conditions that would make the denominators zero.

There are undocumented commands SaveSession and RestoreSession. These are written in top-level Maple in a clean, simple style so they're easy to read, understand, and modify. These save and restore more stuff than just user variables: kernelopts settings, interface settings, debugger settings, enviroment variables, etc. These procedures may not have been updated for several years, and there may be some newer things that are worth saving (like some kernelopts and interface settings) that aren't mentioned in these procedures that you may want to add.

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.