@dharr, sir, The expected results for velocity look like this, means u=1 at x=0 in the range of y=0..1 The actual BCs for u(y) is u(1 + x^2/2) =1  and D(u)(0) = 0, but when we convert these BCs for Psi, it becomes D(psi)(1 + x^2/2) = 1 and D(D(psi))(0) = 0. After converting system in the form of psi by using stream function, the final equation becomes 4th order, and two more Bcs are needed, which are psi(1 + x^2/2) = H, psi(0) = 0.

@dharr sir, here is the final output of my code, and the curves for u(y) started from zero, not 1. 

@dharr thanks I got this point but could you ease resolve issue regarding u(y).

@acer I am solving a system of three coupled ordinary differential equations (ODEs), where the first equation for u(y) contains the unknown pressure gradient term dp/dx​. To eliminate this term and simplify the system, I introduced a stream function such that u(y)=dψ/dy​. This reformulation allows me to rewrite the entire system in terms of ψ(y) and solve for it numerically. Once ψ(y) is determined, compute u(y). However, I observed a problem: according to boundary conditions, the solutions for u(y), θ(y), and ϕ(y) should start from 1, but in my current code, the plots for these variables begin from zero.

The Bcs for u(y) are u(1 + x^2/2) = 1, D(u)(0) = 0, but in terms of psi(y),  is D(psi)(1 + x^2/2) = 1, D(D(psi))(0) = 0.

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