@acer thanks sir for your help kind suggestions....The above file is working well, bt the part 

# ---  Nusselt number for k agaisnt Br --- #
# commented out since it justproduces an error. It's not
# the uiser's main new problem.
(*
BrVals := [0.1, 0.2, 0.3]:
kVals := [0.1, 0.3, 0.5]:  # Modify or extend as needed

# Fixed parameters (except Br and k)
baseParams := [x = 0, We = 0.1, H = 0.65, N[b] = 0.3, N[t] = 0.5]:

# Container for all plots
allPlots := []:

# Loop over different values of k
for j from 1 to numelems(kVals) do
    kVal := kVals[j]:
    NuVals := Array(1..numelems(BrVals)):

  for i from 1 to numelems(BrVals) do
    solNu := dsolve(
        eval(sys, [paramSet[], Br = BrVals[i]]),
        numeric,
        method = bvp[midrich],
        maxmesh = 50000,
        initmesh = 300,
        abserr = 0.1e-7
    ):
    
    # Nu = - dtheta/dy at y=1 (theta is 2nd function)
    NuVals[i] := -evalf(solNu['eval'](1, 2, 1));
end do:

       
    # Convert Nu values to points and plot for this k
    curve := pointplot(
        [seq([BrVals[i], NuVals[i]], i = 1..numelems(BrVals))],
        color = ColorTools:-Color([j/numelems(kVals), 1, 1], colorspace = hsv),
        symbol = solidcircle,
        connect = true,
        legend = cat("k = ", kVal),
        thickness = 2
    ):
    
    allPlots := [op(allPlots), curve]:
end do:

# Combine and show all curves in a single plot
display(
    allPlots,
    labels = ["Br", "Nu"],
    labelfont = ["Times", "italic", 18],
    title = "Nusselt Number vs Brinkman Number for Various k Values",
    titlefont = ["Times", 18],
    axesfont = ["Times", 14],
    size = [800, 600],
    axes = boxed,
    legendstyle = [font = ["Times", 12]]
);
*)
is not producing graphs, and is it possible to save data points for Nu like other quantities? Besides this, the above code is ok for different values of k, but when I run the same code for "We" values, the method is not convergent. is it possible to handle this error?

@acer sir thanks for your help, and it's working well. Besides this, I extended my work and added some extra equations to its solution, and also it seems ok. Actually, I need data points for further analysis, and I even generated data points, but to check those data points, I plotted the graphs against those points, and they are different from the actual plots. It seems there is something wrong. Could you please correct my procedure? Its my last correction please help.

Numeric_help_ac.mw

@acer sir My actual ODE is de1 := diff(u(y), y$2) + We * diff(diff(u(y), y)^2, y) = diff(p(x), x), and the BCs associated with it are u(-1 - x^2/2) = 1, u(1 + x^2/2) = -k. The de1 has diff(p(x), x), and we do not know its value, so we introduced stream function u(y) = diff(psi(y), y) and transformed the whole equation into psi. After that, the equation becomes 4th order and needs 4 BCs, so we convert u(-1 - x^2/2) = 1, u(1 + x^2/2) = -k into psi, i.e., D(psi)(-1 - x^2/2) = 1, D(psi)(1 + x^2/2) = -k, and also we need two more BCs, and for it, 2H = psi(1 + x^2/2) - psi(-1 - x^2/2) from analysis. after comparing it, we have

2H = psi(1 + x^2/2) 

0 = psi(-1 - x^2/2)

My main task is to plot graphs for u(y) against y, and the range for it is -1..1. So, firstly I found the solution of psi(y) and then made a plot for u(y).

 Sir, Could you please help me to fix my issue for the solution of u(y)?

@acer Yes, sir, you are right. The LHS of de5 contains psi(y) and theta(y), which means this part of the equation is dependent on y, and we need to use the solution for psi(y) and theta(y) from the above part of the example in de5. After that, need to plot dp/dx against x for the range -1..1. I am trying to implement this idea on my problem from the following paper.  Please help in this regard.

  https://4spepublications.onlinelibrary.wiley.com/doi/epdf/10.1002/pen.26750

@acer yes sir , you are right; my DE is the function of Psi(y) and theta(y). The relation between psi(y) and u(y) is u(y)= diff(psi(y),y) and I was directly calling u(y) instead of diff(psi(y),y).  Acutally, diff(psi(y),y) means solution of u(y). Thanks for your correction. Beside this, can you please help me to plot diff(p(x),x) agisnt x=-3..3 by fiing other parameters as defined in the file you corrected from Equation (5) by using the solution of psi(y) and theta(y). Also, this equation depends on y, and its value is y = k + (1 - k)*x. Can you please help me regarding the graph of the  Nusselt number by using the solution of theta(y) against any parameter like Br. The procedure for Nusselt is

Nu:= diff(theta, y);
Apply limit on the Nu i.e y=h=k + (1 - k)*x, and then graphs for Nu against Br.

@acer , Sir, Acutally, u(y) is a function of both x and y, and to get results for u(y) against y in the range y=0..1, we suppose to fix a value of x like 0, 0.1, or 0.5. That is why define paramSet early and suppose the value of x=0.  Now to plot the results for diff(p(x),x) against x for range x=0..1 from equation 1 by using the solution of u(y) or from equation 4 by using the results of psi(y).  Two different things: when plotting for u(y), the range for y=0..1 and when needing to plot diff(p(x),x), the range for x=0..1. Acutally, I am following this idea from a paper that has already been published. For reference, I am uploading that paper, and please have a look at equations from 27 to 39. There is some discussion adding on the solution methodology right after eq. 39. Sorry pdf is not supported

@acer , sir can you please have a look on my below comment about the relation between y and x 

@acer Thanks, sir, for your help. The output for u(y) and theta(y) is even suitable but may be I am doing something wrong for diff(p(x),x). I think the problem is from BCs. let me clear rlation between y and x. Actually, in my case, y=h and the value of h is h= k + (1 - k)*x; so from this, y=h= k + (1 - k)*x. For the plot of u(y) and theta(y), we suppose x= 0 in paramSet := [x=0, beta = 0.5, Br=2, m=1.5, A=0.5,H=1.3,k=2.5]: but for diff(p(x),x) we need plot it against x.

my Bcs are 

theta(0) = 1, theta(k + (1 - k)*x) = 0, D(psi)(0) = 1, D(psi)(k + (1 - k)*x) = 0,psi(0) = 0, psi(k + (1 - k)*x) = H}:

In these BCs y=h=k + (1 - k)*x

@acer, sir Eq. (1) is my main equation, which is denoted by de1. I need to solve this equation for u(y), but this equation has diff(p(x), x). So, we define u(y):=diff(psi(y),y) and then eliminate diff(p(x), x) and solve for Psi(y). After that, we solve for u(y) and theta(y).  After getting the solution of u(y), we now need to use the solution of u(y) in Eq. (1) and plot the solution for diff(p(x), x), against x=0..1, which is denoted by dpdx in the maple sheet. There is another way; we can also find the solution for diff(p(x), x) from Eq.(5) by using the solution of Psi(y).  Also, we need to plot p(x) against x by integration or any method on diff(p(x), x) by using the following two conditions: P(0)=0 and P(1)=0. 
Now need to integrate the final solution of p(x) over the interval 0..1 to get numeric values in tabular form. The expected results for diff(p(x), x) and p(x) are .

@mmcdara, @acer , @Kitonum  could you please have look on my above problme?

@mmcdara sir, I replied to your this response, but my reply was automatically adjusted in the above comment.

@mmcdara Yes, sir you are right; the expression PP is dependent only on x. We have to calculate PP_Max from the expression PP by using the second derivative test, and the procedure is defined in the Maple sheet. Then we need to plot PP_Max against k for different values of epsilon. Further, we need to plot PP/PP_Max against X for different values of epsilon. For reference, I uploaded some pictures for both. Our  results should look like that, not exactly the same.

@mmcdara  I checked, and there is no issue in the plot now. Can you please upload the Maple sheet and how to plot PP\PP_max?

Actually, I am trying to find the maximum relation by using the derivative test. My actual expression is a bit complex and difficult to solve, so I use series from means, but this has again a problem. So could you please have a look at my maple sheet and correct the procedure? I take help only but have problems and do not produce results. I explained the procedure in the MPAL sheet, but you do not know whether it is correct or not.

Maximum_help.mw

The results look like this

@Kitonum Thanks sir