@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).

@MaPal93 If I understand you correctly (based on the code in your first Reply), then you consider each column of each matrix as a separate function, so there should be three "curves" for each matrix. Is that correct? If so, does this plot make sense for your purpose?:

X:= <seq(1..op([1,1], M_jk))>:  #op([1,1], ...) is the row count.
plot([seq](<X | M_jk[..,j]>, j= 1..op([1,2], M_jk))); #op([1,2], ...) is the column count.

If that makes sense for your purpose, then do the same for the other matrix, and put the two plots together with dualaxisplot.

If that doesn't make sense for your purpose, try to explain to me why not, and I'll probably have another idea.

@PhearunSeng 

Assuming that 

1. the non-atomic entries (i.e., the inner vectors or matrices) of the outer matrix DOe each contain only 1 row,
   AND
2. the number of atomic entries (the integers in your example) in each row of DOe is the same,

then it can be flattened like this:

#Construct your example:
DOe:= eval(<a,b; b,c; c,a>, [a,b,c]=~ `[]`~([<1|2>, <3|4>, <5|6>]));

#Flatten it:
Matrix(
    [seq](
        (ListTools:-Flatten@[subsindets])(DOe[k], rtable, v-> seq(v)), 
        k= 1..op([1,1], DOe)
    )
);

Maybe now you see that it's far easier to just avoid creating matrices that need to be flattened.
 

@PhearunSeng I don't think that you understood the read-between-the-lines rhetoric of the "flattening" part of my Answer, due to my gentle tone, so I need to reiterate it more vehemently (also reiterating something Tom Leslie said to you):

You have some earlier code, which you haven't shown, which produces the matrix DOFe. That code is a corrupt, filthy piece of crap. That code is what you should ask me to teach you to correct, so that I'd be teaching you to avoid producing crap rather than teaching you how to wipe it up after it's produced.

However, I will show you how to flatten that matrix if you'd attach a worksheet that produces it. Please attach it. The plaintext code that you gave just gives me a error (which is likely not your fault).

Unrelated issue: I see this command in your code:

with(linalg)

This is also filth, crap; "deprecated" is the official word, which sounds vile enough already. You should never use that command.
 

Anecdote: In the mid 1990s, a company that I was doing some computing consulting for purchased a full-package (integrated hardware and software) accounting / inventory / point-of-sale (i.e., cash register) system. The software was written in a ridiculous language named KADOL, which was proprietary to that company. My job was to "maintain" that software, that is, make on-demand and reasonable modifications to it that were peculiar to my client's business (a chain of 11 hardware stores, 7 retail and 4 for building professionals).

All of that KADOL code was crap. It worked, but was very difficult to read and modify. Fortunately, the files of code were each only a few hundred lines long. At the top of one file, where a programmer might ordinarily put their name in a comment, I saw this comment instead:

# Are you kidding? I wouldn't put my name on this crap.

@vs140580 In my formulation, G1 is the smaller graph and G2 the larger. In your followup immediately above, it seems that you may have switched those roles. I'm not sure; your writing is difficult to understand. Can you please check whether in the above G1 always refers to the smaller graph and G2 the larger?

@vs140580 It's okay. Even I didn't know about the importance of the "Text" and "Math" boxes until today. I do believe your sincerity.

@Dkunb Let's suppose that your expression for H[s] is correct:

H[s]:= sum(Norm(p_[n])^2/2/m + U(q_[n]), n= 1..s) 
    + sum(sum(V(q_[i] - q_[j]), j= i..s), i= 1..s)/2;

(I have no reason to doubt that it's correct; I simply say "let's suppose" because my next conclusion depends on it being correct, which I have no means of verifying.) I consider H[s] as a function of positive integer s (function being used here in the mathematical sense, not either of the Maple senses). You want H[N-s]. Via done-in-my-head substitution (s:= N-s) and re-indexing (shifting indices up by s), I get

H[N-s]:= sum(Norm(p_[n-s])^2/2/m + U(q_[n-s]), n= s+1..N)
    + sum(sum(V(q_[i-s] - q_[j-s]), j= s+i..N), i= s+1..N)/2;

Does that seem correct to you? And do you see how it's different (and more than just superficially different!) from what you had?

@nm It can be corrected by changing `if`(e::`=`, lhs-rhs, x->x) to 

(numer@`if`(e::`=`, lhs-rhs, x->x))

And, although it's not necessary for any among your current list of examples, I think it'd be wise to also change evala@Norm to numer@evala@':-Norm'

The final expression (the integral) contains two partial derivatives of "rho", but there's no rho mentioned earlier.