# Integration NDsolve result with a different variable

This is an old question but I am having a new problem. I was able to solve z,r,phi as a function of s but what I really need at the end is to plot z as a function of r. Since dz/dr=Tan[phi], I need to integrate Tan[phi] with respect to r. How should I do this?

With[{Bo = 5},
Eqns = {2 - Bo z[s] - Sin[φ[s]]/r[s] ==
Derivative[1][φ][s],
Derivative[1][z][s] == Sin[φ[s]],
Derivative[1][r][s] == Cos[φ[s]],
z[0] == 0, φ[0] == 0, r[.001] == .001}];

sol = NDSolve[Eqns, {z, r, φ}, {s, 0.01, 3}]
f1 = (Tan[φ[s]] /. sol);
f11 = NIntegrate[f1, {r, 0, 3}];
p1 = Plot[f11, {r, 0.01, 3}, PlotStyle -> Red]


Thanks!

sol = NDSolveValue[Eqns, {z, r, \[CurlyPhi]}, {s, 0.01, 3}]