As discussed in my Reply to Preben, the following code addresses all the issues raised thus far, and does so in a way that's transparent to the end user.

restart:

unprotect(
    LinearAlgebra:-Eigenvectors, LA_EV_orig, LA_EV_inner, LA_EV_sort
):
LA_EV_orig:= eval(LinearAlgebra:-Eigenvectors):

SetOrder:= (x,y)-> x={x,y}[1]
:
LA_EV_sort:= (R::list)->
#R is [output] from LinearAlgebra:-Eigenvectors. The eigenvalues and 
#corresponding vectors are sorted in "set order".
    if R=[] then return 
    elif R=[1]::Vector then 
        local p:= sort(R[1], SetOrder, ':-output'= ':-permutation');
        if nops(R) > 1 and R[2]::Matrix then
            (R[1][p], R[2][.., p], R[3..][])
        else
            (R[1][p], R[2..][])
        fi
    elif R[1]::list then 
        (sort(R[1], SetOrder@op@curry(op,1)~@`[]`), R[2..][])
    else R[]
    fi
:
LA_EV_inner:= proc(AC::seq(specfunc(_Inert_MATRIX))) 
option cache;
    (LA_EV_sort@[LA_EV_orig])(FromInert~([AC])[], _rest)
end proc
:
LinearAlgebra:-Eigenvectors:= overload([
    proc(
        AC::seq(Matrix(square, algebraic)), 
        O::seq({name, name=anything}):= (),
        $
    )
    option overload;
        LA_EV_inner(ToInert~([AC])[], O)
    end proc,

    LA_EV_orig
]):
protect(LinearAlgebra:-Eigenvectors, LA_EV_orig, LA_EV_inner, LA_EV_sort):
    
#Example:
Sy:= Matrix([[0,-I,0],[I,0,-I],[0,I,0]])/sqrt(2):
for k to 5 do (ev||k, EV||k):= LinearAlgebra:-Eigenvectors(Sy) od:
ev1, EV1;

#Prove that they're all equal:
andmap(LinearAlgebra:-Equal, {ev||(2..5)}, ev1), 
    andmap(LinearAlgebra:-Equal, {EV||(2..5)}, EV1);

                       [-1]  [   -1     1     -1    ]
                       [  ]  [                      ]
                       [ 0]  [   (1/2)         (1/2)]
                       [  ], [I 2       0  -I 2     ]
                       [ 1]  [                      ]
                             [   1      1      1    ]
       
                                 true, true

#Request list-form output:
LinearAlgebra:-Eigenvectors(Sy, output= list);
 

The commands (a) evalb(A=B), which compares expressions' identities, and (b) is(A=B), which compares expressions mathematically, possibly under assumptions, are intentionally different, and this is not due solely to efficiency considerations. It wouldn't be "more robust" to change (a) in any way that makes it more like (b). Indeed, it would be far less robust, because (b) must necessarily (see Richardson's Theorem) FAIL is some cases.

The poorly named and mostly redundant command LinearAlgebra:-Equal maps operation (a) over arrays. It's trivial to write a similar command that maps (b) over arrays. There are many possibilities, some shown by the other respondents. My first choice is

Is:= (A::rtable, B::rtable)->
    LinearAlgebra:-Equal(0~(A), 0~(B), _rest) and 
    andmap(is, simplify~(A-~B)=~ 0)
:

This procedure behaves reasonably even if A and B have mismatched sizes or shapes.

Several programming points:

1. The break statement (up until Maple 2021) breaks out of the innermost loop that contains it. Maple 2021 introduced the multi-level break, which is described on pages 1-5 of this document: https://www.maplesoft.com/products/maple/new_features/Maple2021/PDFs/Language.pdf

However, Kitonum has already given the better solution for this particular example: return.

2. The print statement is never an appropriate way for a procedure to return a value. It is only for displaying supplementary information.

3. A procedure should not rely on a particular package having been loaded with with. Nor should a with command ever be used in a procedure. Instead, use a uses clause (example below).

4. The problem of loops nested to an arbitrary depth can usually be solved by using a Cartesian product, like this:

Conremovedges:= proc(
    G::GRAPHLN, 
    Paircrossedges::And(
        {set,list}({set,list}(set)),
        identical({[2,2]}) &under ((nops~@[op])~@{op})
    )
)
uses GT= GraphTheory, It= Iterator;
    for local p in It:-CartesianProduct(Paircrossedges[]) do       
        if 
            (GT:-EdgeConnectivity@GT:-DeleteEdge)(
                G, (local dedges:= {seq}(p)), inplace= false
            ) = 3
        then return dedges
        fi
    od;
    "not found"
end proc
:
Conremovedges(
    GraphTheory:-CompleteGraph(6), 
    [[{1,6},{3,5}],[{2,5},{1,4}],[{2,6},{3,4}]]
);
                    {{1, 6}, {2, 5}, {2, 6}}

​​​​​​

I believe that that is a result of an anti-piracy measure included by Maplesoft in specific version(s) of Maple. It's possible that the measure is triggered by false positives. It has been reported here on MaplePrimes about a dozen times. The solution is to reinstall Maple with correct licensing information. My guesses about all this are based solely on what's been reported here on MaplePrimes, and thus they're highly speculative.

By the way, "Problem with output" is a very poor title for a Question here because it could describe the majority of things that users report.

You need to change self to _self. This is mentioned on page 6 of this PDF: https://www.maplesoft.com/products/maple/new_features/Maple2021/PDFs/Language.pdf

Suppose that your solution is stored in Sol. Then do

NewSol:= eval(Sol, BesselK= ((x,y)-> (2/y)^x/GAMMA(x)));

Note that GAMMA is all capitals.

You might find the command MultiSeries:-asympt useful.

If you "know" f (except for C), then your second relation, f(x^2) + g(x) = sin(x) - x^4, completely determines g (except for C). So, x*D(f)(x) + D(g)(x) is completely determined (the C disappears), and it's not equal to cos(x) - 3*x^3.

Did you perhaps mean, instead of "know", that -x^2/2 + C was just a guess or ansatz for f?

To distribute a denominator over the terms of a numerator, use expand. So, expand(eq_i4).

How about this?

AllEdgesDistance:= proc(G::GRAPHLN, v)
uses GT= GraphTheory;
local 
    VL:= GT:-Vertices(G), V:= table(VL=~ [$1..nops(VL)]),
    D:= GT:-AllPairsDistance(G)[V[v]], e, u
;
    sort([seq]([min(seq(D[V[u]], u= e)), e], e= GT:-Edges(G)))
end proc
:    
P:= GraphTheory:-SpecialGraphs:-PetersenGraph();
AllEdgesDistance(P,1);
[[0, {1, 2}], [0, {1, 5}], [0, {1, 6}], [1, {2, 3}], [1, {2, 9}], [1, {4, 5}],
 [1, {5, 8}], [1, {6, 7}], [1, {6, 10}], [2, {3, 4}], [2, {3, 7}], [2, {4, 10}], 
 [2, {7, 8}], [2, {8, 9}], [2, {9, 10}]]

The counting work can be done---using a single pass through the list---by ListTools:-Collect. 

L:= [A,B,A,C,A,B,D,E,A,F,G,H,H,G,I,P,Q,W,A]:
interface(rtablesize=nops(L)):
MyTable:= <
    <Data | Frequency>,
    (rtable@sort)(
        ListTools:-Collect(L), key= curry(sprintf,"%a")@curry(op,1)
    )
>;

You can then use ExcelTools:-Export to send that to Excel. I don't know about Word.

The set of all vertices at a given distance d (which defaults to 2 in the procedure below) from a given vertex v is returned by this procedure. It avoids the need to scan through all vertices and measure their distance from v. This will be more efficient if the number of vertices is large and d is small.

DistNVertices:= proc(G::GRAPHLN, v, d::nonnegint:= 2)
local V:= {v}, U:= {v};
    to d do 
        V:= `union`(({op}@GraphTheory:-Neighborhood)~(G,V)[]) minus U;
        U:= U union V
    od;
    V
end proc
:

Given any table, you can construct its setwise inverse table with this procedure:

TableInverse:= (T::table)->  (([lhs]~)~@ListTools:-Classify)(rhs, [entries](T, ':-pairs')):

Example of use:

T:= table():  R:= rand(0..5):  for i to 9 do for j to 9 do T[i,j]:= 'R'()$2 od od:

TI:= TableInverse(T):
T[6,7];
                              1, 2

TI[1,2];
                    {[2, 5], [6, 7], [7, 2]}

If the table is one-to-one (aka injective), or you want to ignore multiple indices mapped to the same entry, then a simpler procedure can be used:

TableInverse1to1:= (T::table)-> (table@(rhs=lhs)~@[entries])(T, ':-pairs'):

The time complexity of TableInverse is O(n*) where n is the number of indices, and the * represents unknown but small (< < n) logarithmic factors. The time complexity of TableInverse1to1 is simply O(n).

You can teach diff a recursive formula for the derivative of P(n,x). Once you do that, it'll be able to automatically do higher-order derivatives and/or derivatives of compositions (such as with cos(theta)).

Important: If you want to do this, do not use with(orthopoly). Instead, if and when you want to eliminate P from an expression ex, do eval(ex, P= orthopoly:-P).

`diff/P`:= proc(n,x,z)
    if depends(n,z) then 
        error "derivative of P wrt its 1st argument not defined"
    else
        diff(x,z)*n/(x^2-1)*(x*'P'(n,x) - 'P'(n-1, x))
    fi
end proc
:
#Optional (only affects display): Display P(n,x) with the n as a 
#subscript:
`print/P`:= (n,x)-> nprintf(`#msub(mi("P"),mo("%a"))`, n)(x)
:
dP:= diff(P(n, cos(theta)), theta$2);
simplify(eval(dP, n=2));
simplify(eval(%, P= orthopoly:-P));

I don't know what version of Maple you are using. I am using Maple 2020. So, I don't know if this Answer will work for you. It's useful to make assumptions on x and alpha. If need be, the system can internally derive an appropriate assumption for t based on those. It may be counter-productive to make an assumption on t. 

This works instantaneously for me:

b:= GAMMA(2-alpha)/(1-alpha)/GAMMA(1-alpha):
J:= Int((x-t)^(-alpha)*a*(t-b*ln(1+t/b)), t= 0..x);
value(J) assuming x > 0, alpha > 0, alpha < 1;

If you need to resort to numerical integration, there are appropriate easy ways to do that in Maple. The method that you show is not very good.
 

Here's another way:

numboccur(evalb~(A >~ B), true)