# Evaluation and plot the result of NIntegrate obtained from NDSolve

First I solved coupled ODE using NDSolve. Now I want to integrate the expression (p'''[0] - ((1.0/k)*p''[0]) + (beta*m*m*s[0]) with respect to the parameter x over the interval [0,1]. I need a plot of the integrated expression, say pr as a function of the parameter F. I am trying with the following codes.

a = 0.5;
b = 2.0;
t = 0.0;
Q0 = 1.0;
alp = 0.5;
alpha = 0.5;
beta=1.0;
k = 0.5;
m = 10.0;
h = 1.0 - (a*Cos[Pi*(x - t)]*Cos[Pi*(x - t)]);
eq1 = {s''[y] == (m/alpha)*Sinh[alpha*s[y]]};
bc1 = {s'[0] == 0, s[h] == 1};
eq2 = {p''''[y] - (1/k)*p''[y] + beta*m*m*s'[y] == 0};
bc2 = {p[0] == 0, p''[0] == 0, p[h] == F, p'[h] == 0};
px = ParametricNDSolve[
Join[eq1, eq2, bc1, bc2], {s, p}, {y, 0, h}, {x, F} ];
pr = Evaluate[
NIntegrate[(p'''[y] - ((1.0/k)*p''[y]) + (beta*m*m*s[y])), {x, 0,
1}] /. {y -> 0}]


The graph of Integrated value vs F should look like the following graph. Kindly note that here beta1 is another parameter involved in the problem. But for simplicity, I ignored that parameter. Here we can try for any other parameter that is present in the code itself, for example, the parameter m.

• There are several similar questions posted by you. Instead of giving feedback you prefer to post new questions, no good style! Sep 14, 2023 at 8:44
• I am not getting the desired answer. That's why I am posting. If you are unable to help then please ignore. Thank you. Sep 14, 2023 at 8:55
• If you "don't get the desired answer" you should give a feedback. This interaction makes stackexchange work and honors the effort made by users trying to help you! Sep 14, 2023 at 9:04
• O.K. Thank you. Sep 14, 2023 at 9:10
• How is parameter \[Beta] defined? Sep 14, 2023 at 9:32

ClearAll["Global*"]
a = 0.5;
b = 2.0;
t = 0.0;
Q0 = 1.0;
alp = 0.5;
alpha = 0.5;
beta = 1.0;
k = 0.5;
m = 10.0;

eq1 = {s''[y] == (m/alpha)*Sinh[alpha*s[y]]};
bc1 = {s'[0] == 0, s[h] == 1};
eq2 = {p''''[y] - (1/k)*p''[y] + beta*m*m*s'[y] == 0};
bc2 = {p[0] == 0, p''[0] == 0, p[h] == F, p'[h] == 0};


function ph[...] returns p'''[y] - (1/k)*p'[y] + (beta*m*m*s[y])

ph[h_, F_] :=Block[{sol},
sol = NDSolveValue[{s''[y] == (m/alpha)*Sinh[alpha*s[y]],p''''[y] - (1/k)*p''[y] + beta*m*m*s'[y] == 0,
s'[0] == 0,s[h] == 1, p[0] == 0, p''[0] == 0, p[h] == F, p'[h] == 0},
p'''[y] - (1/k)*p'[y] + (beta*m*m*s[y]), {y, 0, 1} ]
]


integration

xi = Range[0, 1, .02];
hi = Map[1.0 - (a*Cos[Pi*(# - t)]*Cos[Pi*(# - t)]) &, xi];
data = MapThread[{#1, NIntegrate[ph[#2, 1], {y, 0, 1}]} &, {xi, hi}];
ListPlot[data]


addendum( modified with new function ph0[] )

ph0[h_, F_] :=
Block[{sol},
sol = NDSolveValue[{s''[y] == (m/alpha)*Sinh[alpha*s[y]],p''''[y] - (1/k)*p''[y] + beta*m*m*s'[y] == 0, s'[0] == 0,s[h] == 1, p[0] == 0, p''[0] == 0, p[h] == F, p'[h] == 0},
p'''[0] - (1/k)*p'[0] + (beta*m*m*s[0]), {y, 0, 1}]]

Fi = Range[-3, 3];
Fdp=Map[{#, NIntegrate[ph0[1.0 - (a*Cos[Pi*(x - t)]*Cos[Pi*(x - t)]), #],{x, 0,1}]} &, Fi    ]//Quiet
ListPlot[{Fdp, Fdp}, AxesLabel -> {F, dp}, Joined -> {True, False}]


• It seems to be fine. But here you have considered the range of the dependent variable y as [0,1]. Also as I mentioned in my question I want to plot integration of (p'''[0] - ((1.0/k)*p''[0]) + (beta*m*m*s[0]) with respect to the variable x vs the parameter F. However, the desired graph should look different. I have edited the question. Please see. Sep 14, 2023 at 9:22
• Thank you very much. This is the required result. But when I try to run it in my system the code is running for 3-4 hours with lots of errors and I don't get any output. Sep 15, 2023 at 10:31
• See my modified answer. Sorry, don't know we it evaluated yesterday Sep 15, 2023 at 11:05
• Thank you. Now it is working fine. However, I've two questions for my understanding. 1. Why do you consider the range of y, [0,1] instead of [0,h]? 2. What is the reason for writing {Fdp,Fdp} twice in ListPlot[{Fdp, Fdp}, AxesLabel -> {F, dp}, Joined -> {True, False}]? Sep 15, 2023 at 11:50
• Also, is it possible to plot the variation for four different values of the variable m say m=5, 10, 15, 20` in the same figure? Sep 15, 2023 at 12:05