How to nest the output of FindRoot in a table of Derivative evaluations?

Question

I am using ParametericNDSolveValue to numerically solve an ODE system, which I then use to find the definite integral of one of the variables in the solution for different values of a parameter. Here is a simpler (and faster!) version of my code:

eqns={y''[t] + d y[t] == 3 a Sin[y[t]], y[0] == y'[0] == 1};
pfun[din_]:=ParametricNDSolveValue[eqns/.{d-> din}, With[{b=1},b NIntegrate[y[t],{t,0,5}]],{t,0,5},{a}];

which if plotted for three values of "d" looks like this

I am using FindRoot to find the "a" values at which these curves meet the line y=3.

Table[a /. FindRoot[pfun[din][a] == 3.0, {a, 0, 2}], {din, 1, 3, 1}]

What I would like to do next is find the value of the first derivative at each of these points (with the ultimate aim of plotting the derivatives w.r.t. the d parameter values used to make each curve). What I'm stuck on is: do I (1) first generate a list of "a" values using the code above and then somehow plug each element into code like this

Evaluate[Table[Derivative[1][pfun[din]][a], {din, 1, 3, 1}]]

(I'm not clear how to get the derivative function to call a subsequent element in the list for each din value).

or do I (2) embed the FindRoot function somehow into the derivative function? (I don't know how to do this, or do I use Table or Map?)

Answer 1

I could not completely understand your question but I would do the following by changing little bit from your definitions.

eqns = {y''[t] + d y[t] == 3 a Sin[y[t]], y[0] == y'[0] == 1};
pfun[din_] :=ParametricNDSolveValue[eqns /. {d -> din}, 
    With[{b = 1}, b Integrate[y[s], {s, 0, 5}]],{t, 0, 5}, {a}];

Then we will simply use a Module to write function that depending on the value of din will find the derivative at the point where function takes the value of $3$.

fun[din_?NumericQ] := Module[{arg, a},
 arg = a /. FindRoot[pfun[din][a] == 3.0, {a, 0, 2}];
 Evaluate[Derivative[1][pfun[din]][arg]]
 ];
Plot[Evaluate@fun[din], {din, 1, 6}, Frame -> True,FrameLabel -> {d, y'[3]}]

However parabolic type of the derivatives is realistic if you crosscheck the slope of the tangents in the following plot.

Code for the tangent plot:

Clear[tangent];
tangent[din_, eps_] := Module[{arg, der, c, a, cval},
 arg = a /. FindRoot[pfun[din][a] == 3.0, {a, 0, 2}]; 
 der = Evaluate[Derivative[1][pfun[din]][arg]];
 cval = c /. Flatten@NSolve[der *arg + c == pfun[din][arg], c];
 Plot[Evaluate@(der*x + cval), {x, arg - eps, arg + eps},PlotRange -> All, 
  PlotStyle -> Directive[Opacity[.7], Red, Thick, Dashed],
  Epilog -> {Directive[Opacity[.5], Black], PointSize[Large],Point[{arg, 3.}]}
  ]
 ];
tans = Table[tangent[din, .2], {din, 1, 6, 1}];
opt = Epilog -> ( Options[#, Epilog][[1, 2]] & /@ tans);
Show[{Plot[Evaluate@(Table[pfun[i][t], {i, 1, 6, 1}]), {t, 0, 5}, 
Frame -> True]}~Join~tans, opt]

In case this solution does not work for you feel free to leave a comment.