@salim-barzani 

If you want me to help you, first wait until you've received an answer to your question before sending several others.
I'm completely at a loss. 
 

@salim-barzani 

I use Maple 2015 and it's likely that one of the two commands

K := simplify(solve(aux, __K));

den := simplify(denom(K), {eq_8=G(eta)});  #which means: "simplilify wrt the rule {eq_8=G(eta)}

do not perform the same way in your Maple 2021 version.

Can you upload your file in order I can give a look at it?  Maybe I will be able to find a way tio fix this. 
If I'm not you will have to wait for someone who uses maple 2021 to help you.


By the way: Maybe you found it a little unclear why I introduced "constant" K?. Here is an attempt to explain this based on a toy example  Detailed_explanation.mw

@salim-barzani 

It's my pleasure

@Ronan 

I updated the mw file whose end was truncated for some unknown reason (I often observed this while using ImageTools with Maple 2015).
The file containing plot t should be locate in your current directory with name img.jpg.
This is what 

here := cat(currentdir(), "/img.jpg"):
plotsetup(jpeg, plotoutput = here)

do.

Note that 

convert(360, units, 'degree', 'rad') 
                              2 Pi    # Pi and not pi

# More precisely 
convert(360, units, 'degree', 'rad') * Unit('rad')  
                        2 Pi rad

Download Unit.mw

Since I haven't heard from you since I submitted a reply three days ago, I take it you weren't interested. I therefore destroyed its content. 

@AHSAN 

Responses to your reply titled re-check

point [2]
You wrote: 

Now, use the solution of u(y) into Eq. 2 and solve it  for T(y) against y with the help of Bcs given in Eq. 3.
Finally, use the solution of u(y) into Eq. 1 and solve for diff(p(x), x).
 

Equation de1 writes symbolically <A, u(y)> = p'(x) where A is some differential operator and <A, u(y)> denotes its action on u(y). Then p'(x) is obviously equal to <A, u(y)> which is a function of y alone;, let's say F(y).
Thus p(x) = w*F(y) + constant.
 

point [3]
I wrote

Your bc are given at 3 points (-1, 0, 1) which is not supported by Maple.
Either you formulate the problem differently by considering a 1st problem in [-1..0], a second in [0..+1] and interface continuity equations... or you remove on point, for instance 0:

You replied

Ans: I checked the BCs and found correct.
 

Your answer "I checked the BCs and found correct" makes no sense.
Maybe you didn't understand my comment: you cannot solve a 1D boundary value problem by setting boundary conditions at three different points.
Either you remove one of these 3 points, for instance y=0... or you solve two coupled BVPs: one is a BVP in the domain [-1, 0] and the other in the domain [0, +1] and you have to write compatibility (contnuity) conditions at their interface y=0.

___________________________________________________________________________________



Responses to your reply titled look like this

Please do not focus on the plots beat read carefully the conditions I use in my last attached file.
(Note: had you correctly answered my initial comments things could be likely clerer... but you didn't)

I already told you that Br is never defined in you data so I took arbitrarily Br=1, which could easily explain the differences between the T curves.

Concerning the differences between the u-curves : check the data you give in your attached file to those used to produce these last results.
As I simply use dsolve there is no reason ti think that the u(y) curves I provide I wrong

___________________________________________________________________________________

Responses to your reply titled  Solution verifiaction

[1]
When you derive equation 5, whose rhs = diff(p(y), y), wrt y,  you det a rhs  equal to 0: doesn't this mean that diff(p(y), y) = 0 and then that p(y)=a+b*y ?

[2]
Note that you write lhs(eq 5) = diff(p(x), x) which I thought fristly (and maybe wrongly) was a typo. The reason us that the lhs of eq 5 is a function of y alone while its rhs is a function of x alone.
To get de4 you differentiate de3 wrt y , which gives rhs(de4) = 0 .

dsolveing  de4 then givesautomatically a constant value for diff(p(x), x).
Look  CODE REMOVED BY THE AUTHOR  see here  for an explanation on a toy example.

So it is completely impossible to retrieve your second plot: it is not something that can be deduced from your equations nor from functions psi(y) and T(y)

___________________________________________________________________________________

Here is a revised version of my work   CODE REMOVED BY THE AUTHOR  see here  

You are probably going to say that the results this file contains are not the same than those you presented before. Nevertheless I advice you to compare the equations you wrote to those where "your" results come from before posting a reply.
For now on you presented us a problem which is not correctly (well?) posed.

Clicking on the link returns "page not found". Likely because of non latin characters in its name.

@AHSAN 

        CODE REMOVED BY THE AUTHOR  see here

ALL CODE SNIPPETS DELETED BY THE AUTHOR  see here

[1]
What does diff(p(x), x)  mean in equation 5?
Wouldn't it be diff(p(y), y)  instead? But if it is so you cannot solve equation 5 numerically as p(y) is undefined.
  

[2]
You have 2 odes (equations 5 and 6) that psi(y) must verify... you must decide the one you want to keep unless they both have the same solution, whixh is highly unlikely given thefirst one depends
on p.

[3]
Your bc are given at 3 points (-1, 0, 1) which is not supported by Maple.
Either you formulate the problem differently by considering a 1st problem in [-1..0], a second in [0..+1] and interface continuity equations... or you remove on point, for instance 0


[4]
I psi(y) verifies equation 6 and bcs 7 you can plot u(y) this way

        SNIPPET DELETED BY THE AUTHOR  see here

And if you want to plot u(y) for different values of k:

​​​​​​​[5]I psi(y) verifies equation 6 and bcs 7 you can plot the rhs of de1 this way
  
​​​​​​​

If it is p(y) that you want to plot, its expression is given at the end of point [3] for k=0.


Reformulate your problem correctly and come to us when it's done

@C_R 

I found this inadvertedly and I use it a lot, but you are right, it is likely an undocumented feature.
I guess the option is accepted because its syntax is the same as PLOT's

Compare these:

plot(1, x=0..1, color=blue); # documented
op(%);  #option 1

plot(1, x=0..1, color=COLOR(RGB, 0, 0, 1)); # undocumented?
op(%); #same structure than option 1

plot(1, x=0..1, color=ColorTools:-Color([0, 0, 1]));  # documented in plot/color
op(%); #same structure than option 1

@vv

It seems that the orthogonal polynomials which ease the interation the most are Chebyshev U.
If the solution of order n is S_n(x) = C₀*U(0, x) + ... + C_n*U(n, x), then the solution of order (n+1) should be 
S_n(x) = S_n(x) + C_n+1*U(n+1, x).

More of this the results are significantly better when numeric integration is done instead of exact integration(???)

Galerkin_solution.mw

@vv 

In fact, I had something else in mind that I didn't see through to the end.

Using othogonal polynomials enables building successive approximations but just computing a single new term for each refinement.
So I began coding this with Hermite polynomials in mind, while I realized that the singulatity of

phi(y)/(x-y+2)^2

will cause problems to compute 

Int(H(m, y) / (x-y+2)^2 *  exp(-y^2), y = -infinity..+infinity)

So I resricted th integration range in order to avoid this singularity but at the same time I unfortunately kept the Hermite poynomials.

I might as well have written

Galerkin := proc(EQ, M, R, d)
  local G, S, N, r, m, n, C:
  G := Matrix((M+1)$2):
  S := Vector(M+1):

  r := evalf(R):

  Digits := d:
  for m from 0 to M do
    N := int((x^m)^2, x=r);
    for n from 0 to M do
      G[m+1, n+1] := int( (value(eval(lhs(EQ), phi = (u -> u^n))) assuming x > 0) * x^m, x=r) / N
    end do:
    S[m+1] := int(rhs(EQ)*x^m, x=r) / N;
  end do:
  G := evalf(G):
  S := evalf(S):
  C := LinearAlgebra:-LinearSolve(G, S):

  return add(C[i+1]*x^i, i=0..M)
end proc:

which gives (obviously) the same final result (R=0..1 for instance)

Thank for pointing out this error.

@Carl Love 

but it's too top-notch programming for me to be able to appreciate all its subtleties.

@Carl Love