# Integration inside and outside of NDSolve

EQU = {Y'[s] == -2 X[s] /Sin[ 2 Y[s]] Z[s],
X'[s] == Sin[X[s]] Cos[Y[s]]/3, Z'[s] == Cos[Y[s]] Cos[X[s]],
X == 0, Y == .3, Z == 1, V1'[s] == Z[s]^2 Cos[Y[s]],
V1 == 0};
NDSolve[EQU, {X, Y, Z, V1}, {s, 0, 2}];
{x[u_], y[u_], z[u_], v1[u_]} = {X[u], Y[u], Z[u], V1[u]} /. First[%];
Plot[{x[s], y[s], z[s], v1[s]}, {s, 0, 2}, PlotLabel -> "Z X Y V1"]
v2[s_] = NIntegrate[z[s]^2 Cos[y[s]], {s, 0, 2}];
Plot[{x[s], y[s], z[s], v2[s]}, {s, 0, 2}, PlotLabel -> "Z X Y v2"]


EDIT1:

V1 is calculated inside NDSolve correctly but the same function (forgotten in the first instance) is calculated with a different name V2 outside NDSolve separately. I used earlier obtained function definitions for constituents as input in NIntegrate. But it does not work. There is something I miss in V2 ( where V1 calculated correctly from initial point of differential equation solution, but fails for V2.) Please help to find it out.

EDIT2

r1 = 3.25; ro = 1; lmax = 5;
EQN := {ph1'[
l] == -Cos[ph1[l]]/
rr1[l]*(2 rr1[l]^2 - 3 ro^2)/(rr1[l]^2 - ro^2),
rr1'[l] == Sin[ph1[l]], z1'[l] == Cos[ph1[l]], ph1 == 0,
zz1 == 0, rr1 == r1, vol1[l] == Pi rr1[l]^2 Cos[ph1[l]],
vol1 == 0};
NDSolve[EQN, {ph1, rr1, z1, vol1}, {l, 0, lmax}];
{r[u_], ph[u_], z[u_],
vol[u_]} := {rr1[u], ph1[u], z1[u], vol1[u]} /. First[%];
P1 = ParametricPlot[{z[l], r[l]}, {l, 0, lmax},
GridLines -> Automatic, AspectRatio -> Automatic]
ParametricPlot[{z[l], ph[l]}, {l, 0, lmax}, AspectRatio -> 1]
volume[smax_] := NIntegrate[Pi r[l]^2 Cos[ph[l]], {l, 0, smax}]
Plot[{ph[l], volume[l]}, {l, 0, lmax}, AspectRatio -> 1]


The above does not work, although earlier it was ok.

EDIT3

Error messages before correct results suggested by MMM remedy for EDIT2: • Could you be MUCH clearer about what exactly your problem is? Feb 7 '17 at 19:02
• Hope it is more clear now. Feb 7 '17 at 20:31

One way to do it is

{Y[s_], Z[s_], V1[s_]} = {Y[s], Z[s], V1[s]} /.First@NDSolve[EQU, {X, Y, Z, V1}, {s, 0, 2}];
V2[y_] := NIntegrate[Z[s]^2 Cos[Y[s]], {s, 0, y}];
Plot[{V1[s], V2[s]}, {s, 0, 2}, PlotRange -> All,
PlotStyle -> {Green, Directive[Red, Dashed]}, Frame -> True] Please avoid using single capital alphabets.

Response to Edit 2

zz1==0 should be z1==0. vol1[l] should be vol1'[l].

r1 = 3.25; ro = 1; lmax = 5;
EQN = {ph1'[l] == -Cos[ph1[l]]/ rr1[l]*(2 rr1[l]^2 - 3 ro^2)/(rr1[l]^2 - ro^2),
rr1'[l] == Sin[ph1[l]], z1'[l] == Cos[ph1[l]], ph1 == 0,
z1 == 0, rr1 == r1, vol1'[l] == Pi rr1[l]^2 Cos[ph1[l]],vol1 == 0};
{rr1[l_], ph1[l_], z1[l_],vol1[l_]} = {rr1[l], ph1[l], z1[l], vol1[l]} /.
First[NDSolve[EQN, {ph1, rr1, z1, vol1}, {l, 0, lmax}]];
volme[x_] := NIntegrate[Pi rr1[l]^2 Cos[ph1[l]], {l, 0, x}]
Plot[{vol1[l], volme[l]}, {l, 0, 5}, AspectRatio -> 1] • It works but with error messages.. Feb 9 '17 at 20:58
• nice answer, note that you can extract the solution functions also like this: {X, Y, Z, V1} = NDSolveValue[EQU, {X, Y, Z, V1}, {s, 0, 2}]. Feb 9 '17 at 21:06