@Mariusz Iwaniuk 

I have detected where the problem lies. The Maple 2020 you used gave a different solution to the same system of equations, which is the correct one, from what Maple 2016 that I am using dished out. This is really serious and disturbing! I hope those in charge can do something about it.

Thanks once again Iwaniuk, at least I have the correct solution now! :)

M2016_vs_M2020.mw

@Mariusz Iwaniuk

It seems to be Maple 2016 issue, I got same problem after executing your file.  Thanks for the effort

@Mariusz Iwaniuk Thanks. Simplication is not the issue, if you look at the last component in the results of the pdetest, the functional derivatives are different and thus cant cancel out. Though I still try your suggestion but got same problem as posted earlier.

solution_bug_1.mw

@Carl Love You are absolutely correct. Just want to make sure it is my Maple version problem. Thanks

@acer Ofcourse it does, and that is what is used for Maths expressions. Want I meant is I just want the substitution ONLY without evaluating the expression as shown in my second attachment. For instance, I want (A_x)_t and not A_xt. The later is what the eval command does.

@acer Thanks, but I actually need the substitution.

As for your question, I was just playing around with different Maple commands when I got the error, and it actually works for an equation I tried earlier as shown in the attached file.

Jusy_trying_stuffs.mw

@Carl Love So far there is no any issue with the equations I have tested, the code is handling all the polynomial nonlinearities effectively. Nevertheless I will take note on the change of the  frontend argument. As for u[x]*sin(u[x]) and its relatives, luckily they are not on the menu for now.

@Carl Love It is ok now. I dully appreciate your resolute (for the past 4 days) to seeing my problem to desired conclusion. Thanks to @acer too for resolving the version problem. Remain blessed!

@Carl Love You nailed it, It works perfectly! The coefficients are exactly correct  but there is an error displayed just before the results in Eq. (4). And also, I don't know why the (----)u_x (same happens depending on the term put in Cfs[...]).

Thanks infinitely many times.

Carl_code_all.mw

@Carl Love  Both codes are working perfectly except for coefficients of u_xu_y,..etc. Thanks for all your efforts.

@Carl Love 

Apologies for my english. What I mean is coefficients of product like:

u_xu_y, u_xxu_xyz, u²_ytu_xyztu_zz..etc.

How can I adapt Carl's code to general settings?  For instance in Eq. (3) of the attached file, how to automatically extract all the coefficients of:

1. u_ijkl,v_ijkl (or their multiples) while the rest are considered constant?

2. each derivatives of the dependent variables and their multiples up to the last order.

I am reposting this question again following @acer suggestion although @Carl Love has replied to my queries, I was aboout to reply to that but found out the thread has been deleted. Carl's answer was to adapt

Jxy := indets(wave, `~`[myindex](`~`[op](0, Jets), identical(`~`[op](Jets)[])$2))

the way one want. Yes I know that because it was stated clearly in the original codes. But I could only achieve that by listing all the derivatives manually as Carl suggested in the deleted thread which may be untractable for large output which is why I mentioned "...how to automatically..." in my question. As for the linearity of the PDEs for code extraction, that is not necessary because what is sought is to create a linear system of PDEs from a nonlinear one by seperating the coefficents of the derivatives of the unknown (and their multiples)  just as been done in obtaining Determining Equations of DEs.

Thanks.

 Adapt_Carl_code.mw

@acer You are right, will do that right now

@Carl Love Yes I know that is where you difined it, you stated it clearly in your code. I tried to manipulate that fancy code but failed, and writing

Jxy:= {u[x,x], u[x,y], u[y,y], v[x,x], v[x,y], v[y,y].....}

may be cumbersome for large equation such as the one in the attached file above which is why I said:

"... how to automatically extract all the coefficients of:"

As you rightly stated, the extractions a meant for linear DEs.

Thanks in advance!

@Carl Love , @acer . Thanks to both of you for coming to my rescue. I have been scratching my head to detect why Carl's code was not working for me. The revision works perfectly.