@janhardo sir, attention is to view nuemric values in tabular form at one point of y such as 0.1 for each k=0.1,0.3,0.5,0.7,0.9 and also looking for to generate 1001 data points of skin friction in the range of y=-1...1 for each k values.

@janhardo Thanks, sir, for your help, but in my system, there is an error.

@janhardo, Thank you, sir; it's working perfectly. I had also submitted another question regarding how to export data points of a different quantity (from the same system) to an Excel sheet. However, it seems that my question was deleted by the system, possibly because it was marked as a duplicate.

I need further assistance on how to generate tables for the skin friction parameter and correctly save the data points of another quantity into Excel. Although the line graphs are generating correctly, there is an error when exporting the data to Excel. Also, in addition to line graphs, is it possible to create 3D graphs for this quantity? Please have a look at my Maple worksheet.

Data_points_help.mw

@janhardo Good day, sir, and thank you very much for your help. In the first part of the file, the tables and graphs of the Nusselt number might have been skipped by you, and they were also showing errors in my file, which I uploaded. If you don't mind, could you please help me with that as well?

@janhardo Dear sir, Thank you very much for your detailed responses and kind assistance. I sincerely appreciate the way you have explained everything thoroughly. Following your guidance, I attempted to run a different system; however, I am currently facing issues with the boundary conditions, and the results are not being generated properly. Could you please assist me with this?

Additionally, is it possible to export 1001 data points of the Nusselt number to an Excel sheet at a desired place in the computer? Specifically, I want the first column to represent the Brinkman number (Br) ranging from 0 to 1, and the second column to contain the corresponding Nusselt number (Nu) values for each fixed value of kkk, such as k=0.1, 0.3,…, 0.9.

Nusselt_Number_help.mw

@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