@Preben Alsholm Thank you. That second example  is more concerning to me because there's no way to know which commands one needs to apply forget to. The unwithing routine for Physics should have the responsibility of taking care of the forgetting.

The example that led to my discovery of this anomaly was a Question here on MaplePrimes that I was attempting to answer yesterday. The user had done with(Physics) just to use Physics:-Setup(...). Then we tried to remove some higher-degree terms from a complicated PDE containing some derivatives of conjugates of the unkown functions. After many attempts, I noticed that the result was completely different depending on whether with(Physics) had been used. Naturally, the next thing I tried was replacing diff with :-diff, to no avail.

By the way, I have no opinion about which form of the derivative of conjugate is better. I just think that there should be a way to switch from one to the other without rigamarole.

@acer Thank you. I guess that there's a bug in forget when the syntax forget(.., whatever) is used. 

The procedure `simpl/conjugate` is on the list of 130 procedures that I just posted above.

@vv For what it's worth, here is a list of 130 top-level-named procedures and appliable modules that reference either the Physics package or global `Physics/Initialized`. (Yes, I realize that many of these are members of older table-based packages and are thus in some sense not top-level names.)

restart:
kernelopts(opaquemodules= false):
interface(warnlevel= 0):
Ph:= CodeTools:-Usage(
    select(
        P-> local pn;            
            try        
                (pn:= ``||(P))::callable
                and indets(
                    ToInert(eval(`if`(pn::procedure, pn, pn:-ModuleApply))),
                    identical("Physics", "Physics/Initialized")
                ) <> {}
            catch: false
            end try,        
        index~(op~(1, {LibraryTools:-ShowContents()[]}), ..-3)
    )
): 
nops(Ph);
(p-> printf("%A\n", p))~(sort([Ph[]])):  
memory used=0.52GiB, alloc change=211.81MiB, 
cpu time=11.24s, real time=10.50s, gc time=2.25s

                              130

ODEtools/Xcommutator
ODEtools/Xgauge
ODEtools/ch
ODEtools/collect
ODEtools/equinv
ODEtools/firint
ODEtools/firint/PDE
ODEtools/first_integral
ODEtools/get_diff_vars_in_expression
ODEtools/info
ODEtools/isolate
ODEtools/odepde
ODEtools/odsolve
ODEtools/redode
ODEtools/remove/hard
ODEtools/subs/do
ODEtools/y(x)_to_y
PD/convert/compseq
PD/diff
PDEAdvisor/high_order/Fourier
PDEAdvisor/high_order/Hankel
PDEAdvisor/high_order/Laplace
PDEAdvisor/high_order/Mellin
PDEtools/U_to_u
PDEtools/difforder/do
PDEtools/isolate
PDEtools/pdesys
PDEtools/useD
PDEtools/useD/diff_to_D
PDEtools/useD/onediff
assuming
assuming/names
combine/cmbtms
combine/exp
combine/linear
combine/linear/do
combine/linear/sort
combine/power
combine/wrt
convert/D
convert/DESol
convert/eval/from_diff
dchange/diff/tr+itr
dchange/funcs
dsolve
dsolve/IC
dsolve/IC/easy
dsolve/IC/get_x_val_and_diff_ord
dsolve/IC/implicit
dsolve/PDE
dsolve/func_in_args
dsolve/general_and_particular
eval/Diff
eval/diff
evalc/abs
evalc/evalc
expand/conjugate
expand/exp/+
expand/ln
expand/power
index/PhysicsTensor
intat
intfactor/formal
is/internal
is/internal/eqns_inds
is/preproc
is/solve
latex
liesymm/convert/D
odeadv/MRL/ch
odeadv/sym_pat/[F(x)+G(y),0]
odetest
odetest/ODESolStruc
odetest/linear
odetest/series
odsolve/Riccati/symmetric
odsolve/high_order/exp_reduce/integrate
odsolve/parametric/Lie
odsolve/solvable/do
odsolve/solvable/linear_homogeneous
odsolve/sort_ODESolStruc/DESol_reduction
pdetest
pdsolve
pdsolve/BC/2nd_order/Series
pdsolve/LagrangeCharpit
pdsolve/ToMissingDependentVariable
pdsolve/checksol
pdsolve/fixedfunc
pdsolve/iPDE
pdsolve/info
pdsolve/resinco
pdsolve/sep/casesplit/do
pdsolve/strip/solve/dsys
pdsolve/sys/hint2
pdsolve/sys/info
pdsolve/sys/sequential
pdsolve/sys/sort
pdsolve/sys/sort/lex
print/Diff
print/Int
print/diff
print/int
property/ObjProperty
redode/many
sdsolve/odeset
sdsolve/ok
sdsolve/process_answer
sdsolve/sequential/sys_and_funcs
sdsolve/sysODE/islinear
signum
simpl/abs
simpl/conjugate
simplify/commonpow
simplify/conjugate
simplify/table
symgen/formal
tools/RewriteRoutinesToUsePhysicsCommands
tools/Typeset
tools/UsePhysics
tools/UsePhysicsOperators
tools/atomize_names
tools/degree
tools/df_dt
tools/map
tools/operate_on_sequence
tools/type
type/PDEtools/known_function
type/PhysicsUnitVectors
type/PhysicsVectors
value/Intat/do

@Ahmed111 In Maple 2022, the result depends on whether or not you use with(Physics). If you don't, the result is 0. So, in your Maple 18, it may be 0 even if you do use with(Physics). I note (see the Question that I recently posted) that Physics:-diff and regular diff treat conjugate differently, and that once either with(Physics) or Physics:-diff is used, it's impossible to get back the regular diff behavior without a restart. 

@Ahmed111 How do you want the derivatives of the epsilons to be handled in the linearization? For example, I could remove all terms that contain an epsilon derivative (of any order) raised to power 2 or higher.

@Rouben Rostamian  Your solution linearizes with respect to x and t, but doesn't that remove from the expression more content than the OP intended to be removed? Can we limit the linearization to the epsilons?

@MaPal93 What you are missing is that you have an extra restart command; thus, matrices K1 and K2 have been erased from memory.

I totally agree with mmcdara about Document Mode and 2D Input. Each is independently an unusable abomination, so the hideousness of the combination is essentially the product of their individual hideousnesses.

@ecterrab Why should there be a way to avoid using the name quotes that is specific to the "prime" character? It's too ad hoc. What about other non-alphanumeric characters? The tremendous number of these ad hoc solutions added in recent years will lead to Maple being unmaintainable as a software product.

By using the name quotes, the ' character can be used in names and as a differentiation operator within a single expression. 

The ASCII (keyboard) character "aposthrope" (') is not even the correct "textbook" character for a "prime" symbol. Observe the difference in these two variations of "f-prime":

`f'`, `#mrow(mi("f"), mo("′"));`;

The one on the left is hideous---I'd never let my work be seen in public like that (with characters overlapping). The one on the right is much closer in appearance to a properly typeset textbook.        

@Ahmed111 Sorry, I missed seeing the derivatives of the epsilon functions (due to the original having used alias to suppress the (x,t) part). Yes, that certainly changes the approach to linearization required. It's a little bit more complicated, but not much more. I'll fix it.

@Ahmed111 After doing mtaylor, do

Eq3:= subs([epsilon[1]= epsilon[1](x,t), epsilon[2]= epsilon[2](x,t)], Eq2):
alias(epsilon[1]= epsilon[1](x,t), epsilon[2]= epsilon[2](x,t)):
Eq3;
 

Side-by-side plots are made by passing a row vector to plots:-display, with each element of the vector being a plot.

#This code assumes that the 2 matrices have the same dimensions.
X:= <seq(1..op([1,1], A_jk))>:    
plots:-display(
    Vector[row](
        op([1,2], A_jk),
        j-> plots:-dualaxisplot(
            plot~([<X | A_jk[.., j]>, <X | A_ki[.., j]>], color=~ [blue, red])[],
            symbol= [solidcircle, cross, asterisk][j]
        )
    ),
    axes= boxed, thickness= 0.2, style= pointline, symbolsize= 8
);

I don't know how to copy-and-paste side-by-side plots into MaplePrimes, but you'll see them in your worksheet. 

By the way, splitting a matrix into column vectors is trivial, and it's done above via the [.., j] syntax.

I get this

P:= proc(M::Matrix) 
local j, X:= <seq(1..op([1,1], M))>;
     [seq](<X | M[.., j]>, j= 1..op([1,2], M))
end proc:
    
plots:-dualaxisplot(
    plot~(
        P~([A_jk, A_ki]), 
        color=~ [blue, red], symbol= [solidcircle, cross, asterisk]
    )[],
    axes= boxed, thickness= 0.2, style= pointline, symbolsize= 8
);

And you can add axis labelling, title, legend, etc. The problem with the above plot is that the magnitude of the first-column data overwhelms the viewing of the other columns. So, how about 3 side-by-side dualaxisplots, one for each column?

Okay, I retract the first sentence of my Answer. Now I think that you do understand what dualaxisplot is supposed to do. What you didn't understand was the correct matrix indexing needed to achieve it.

@MaPal93 I don't see the point in switching the matrix from A_jk to M_jk with some intermediate structure named R_jk. The immediate problem is that R_jk isn't defined in your most-recent worksheet, which makes M_jk a matrix full of symbolic R_jk nonsense. Just use A_jk:

X:= <seq(1..op([1,1], A_jk))>:
plot([seq](<X | A_jk[..,j]>, j= 1..op([1,2], A_jk)));

@mmcdara Actually, it was unclear to me whether dualaxisplot would ultimately be correct for the OP's purpose, because I didn't know what the purpose actually was. Whatever the purpose, the individual plots must make sense before merging them with dualaxisplot (as I'm sure you're aware).