# How to take the derivative of a NDSolve solution?

I have solved and plotted a system of differential equations, shown below, using Mathematica's NDSolve function. Now I would like to take the derivative of one the resulting interpolating functions and evaluate it at a point in the domain.

Pbub = 1.00002
kl = 2.09015;
ccl20 = 0.00206957
cc2h4eq = 0.000301892
Dc2h4 = 1.48*10^(-5);
Dcl2 = 2.90*10^(-9);
t = 183.0 + 273.15;
k1 = 11493.0*Exp[-2156.58/t];
k2 = 8.517*10^9*Exp[-7282.21/t];
\[Delta] = Dc2h4/kl;
eqn1 = Dc2h4*aa''[x] == k1*aa[x]*b[x] + k2*aa[x]*b[x]^2;
eqn2 = Dcl2*b''[x] == k1*aa[x]*b[x] + 2*k2*aa[x]*b[x]^2;
eqn3 = Dcl2*c''[x] == -k1*aa[x]*b[x];
eqn4 = Dcl2*d''[x] == -2*k2*aa[x]*b[x]^2;
eqn5 = Dcl2*ee''[x] == -2*k2*aa[x]*b[x]^2;
bc1 = aa == cc2h4eq;
bc2 = aa[\[Delta]] == 0;
bc3 = b' == 0;
bc4 = b[\[Delta]] == ccl20;
bc5 = c' == 0;
bc6 = c[\[Delta]] == 1.253/98.95;
bc7 = d == 0;
bc8 = d[\[Delta]] == 0.0001;
bc9 = ee == 0;
bc10 = ee[\[Delta]] == 0.0001;

soln = NDSolve[{eqn1, eqn2, eqn3, eqn4, eqn5, bc1, bc2, bc3, bc4, bc5,
bc6, bc7, bc8, bc9, bc10}, {aa[x], b[x], c[x], d[x], ee[x]}, {x,
0, 4*\[Delta]}]


Specifically, I would like to determine: aa'. I've tried various methods using assignment operators, evaluate, etc ... However, none of them have worked thus far. How should I proceed?

• Does the above code run without errors for you? I get error NDSolve::ndsv: Cannot find starting value for the variable aa^\[Prime]. You should make sure code posted runs without errors. You do not have Pbub defined. Nov 16, 2017 at 6:22
• I left out code that defined a few of the required parameters. The post has now been edited to include them.
– Tunk
Nov 16, 2017 at 6:28
• You can apply the differential operator D (see documentation) also to InterpolatingFunctions. Nov 16, 2017 at 6:31
• @HenrikSchumacher I attempted to use D, however, I could not get it to work.
– Tunk
Nov 16, 2017 at 6:34
• Include aa' in your list of functions to be solved for. Nov 16, 2017 at 6:46

First, better to use aa and not aa[x]

soln=NDSolve[{eqn1,eqn2,eqn3,eqn4,eqn5,bc1,bc2,bc3,bc4,bc5,bc6,bc7,bc8,bc9,bc10},
{aa,b,c,d,ee},{x,0,4*\[Delta]}]


Now

Plot[Evaluate[aa[x]/.soln],{x,0,1}] The derivative is D[aa/.soln] And

D[aa[x]/.soln,x]/.x->0 • Thanks a ton! What is the problem with regards to using aa[x] versus aa?
– Tunk
Nov 16, 2017 at 6:39
• @Jay there is really not a problem per say, but with NDSolve it is much easier to use aa instead of aa[x] in the call (vs. with DSolve where aa[x] would be the better choice there. It just makes post processing easier of the interpolation functions, and most of the examples in help are also this way. Nov 16, 2017 at 6:40