Return to Answer

As per comments below, FindFormula only works with one variable. So the following uses FindFit and applies data of myIntFun[r,1,2,$\phi$] to: $$ a+br+cr^2+d\phi+e\phi^2+fr\phi $$ in the range of $0.1\leq r\leq 0.9$ and $0.1\leq \phi\leq 0.9$ and compares the fitted formula to a ListPointPlot3D of the data points:

As per comments below, FindFormula only works with one variable. So the following uses FindFit and applies data of myIntFun[r,1,2,$\phi$] to: $$ a+b\phi+c\phi^2+dr+er^2+fr\phi $$ in the range of $0.1\leq r\leq 0.9$ and $0.1\leq \phi\leq 0.9$ and compares the fitted formula to a ListPointPlot3D of the data points:

Lots of functions to fit numerical data (I use FindFormula below). Below I start with setting up the integral function with integration limits of $1\leq z\leq 10$ for starters and vary just $\phi$ generating a table of {$\phi$,integralFunction[r,p,d,$\phi$]}. Then use FindFormula to fit this data.

**** See update as per comments below ****

Lots of functions to fit numerical data (I use FindFormula below). Below I start with setting up the integral function with integration limits of $1\leq z\leq 10$ for starters and vary just $\phi$ generating a table of {$\phi$,integralFunction[r,p,d,$\phi$]}. Then use FindFormula to fit this data.

Update:

As per comments below, FindFormula only works with one variable. So the following uses FindFit and applies data of myIntFun[r,1,2,$\phi$] to: $$ a+br+cr^2+d\phi+e\phi^2+fr\phi $$ in the range of $0.1\leq r\leq 0.9$ and $0.1\leq \phi\leq 0.9$ and compares the fitted formula to a ListPointPlot3D of the data points:

phiTable = Table[ {phi, r, myIntFun[r, 1, 2, phi]}, {phi, 0.1, 0.9, 0.05}, {r, 0.1, 0.9, 0.05}]; myParms = FindFit[Flatten[phiTable, 1], a + b x + c x^2 + d y + e y^2 + f x y, {a, b, c, d, e, f}, {x, y}]; myFit[x_, y_] = (a + b x + c x^2 + d y + e y^2 + f x y) /. myParms; my2DPlot = Plot3D[myFit[x, y], {x, 0, 1}, {y, 0, 1}, PlotStyle -> {Opacity[0.2], Blue}]; lp = ListPointPlot3D[phiTable, BoxRatios -> {1, 1, 1}]; Show[{lp, my2DPlot}] 

Lots of functions to fit numerical data (I use FindFormula below). Below I start with setting up the integral function with integration limits of $1\leq z\leq 10$ for starters and vary just $\phi$ generating a table of {$\phi$,integralFunction[r,p,d,$\phi$]}. Then use FindFormula to fit this data.

Then I would try to vary two variables like $\phi$ and $r$, generate a table of {$\phi$,r,integralFunction[r,p,d,$\phi$]} and then try to fit this 3D data. Then try to generate a table of another variable and so on. Kinda brute force but it's something you may wish to work with.

Here is a start just varying $\phi$ from 0.01 to 0.99. The red fitted function fits the points nicely.

(* define functions *) wFun[z_, r_, p_, d_] := ( r (d + (r + z) Sech[(z + p)/d]^2 - d Tanh[(z + p)/d]) (z + r Tanh[(z + p)/d]))/(d (r + z)^3 Sech[p/d]^2); qFun[z_, r_, \[Phi]_] := -2 Exp[( z \[Phi] (6 r^2 (\[Phi] - 1)^2 - 3 r z (\[Phi]^2 + \[Phi] - 2) + 2 z^2 (\[Phi]^2 + \[Phi] + 1)))/(2 r^2 (\[Phi] - 1)^3)]; (* define integral function *) myIntFun[r_?NumericQ, p_?NumericQ, d_?NumericQ, \[Phi]_?NumericQ] := NIntegrate[wFun[z, r, p, d] qFun[z, r, \[Phi]], {z, 0, 10}]; (* for now, create table varying just phi from 0.01 to 0.99 and fit a \ formula to this 1D data. *) phiTable = Table[ {phi, myIntFun[1, 1, 2, phi]}, {phi, 0.01, 0.99, 0.01}]; (* find a formula for data *) theF[x_] = FindFormula[phiTable, x] (* superimpose ListPlot of points with fitted function *) lp = ListPlot[phiTable, PlotStyle -> {Black, PointSize[0.01]}] p1 = Plot[theF[x], {x, 0.01, 0.99}, PlotStyle -> Red] Show[{lp, p1}] 

Lots of functions to fit numerical data (I use FindFormula below). Below I start with setting up the integral function with integration limits of $1\leq z\leq 10$ for starters and vary just $\phi$ generating a table of {$\phi$,integralFunction[r,p,d,$\phi$]}. Then use FindFormula to fit this data.

Then I would try to vary two variables like $\phi$ and $r$, generate a table of {$\phi$,r,integralFunction[r,p,d,$\phi$]} and then try to fit this 3D data. Then try to generate a table of another variable and so on and also increase the integration limits. Kinda brute force but it's something you may wish to work with.

Here is a start just varying $\phi$ from 0.01 to 0.99. The red fitted function fits the points nicely.

(* define functions *) wFun[z_, r_, p_, d_] := ( r (d + (r + z) Sech[(z + p)/d]^2 - d Tanh[(z + p)/d]) (z + r Tanh[(z + p)/d]))/(d (r + z)^3 Sech[p/d]^2); qFun[z_, r_, \[Phi]_] := -2 Exp[( z \[Phi] (6 r^2 (\[Phi] - 1)^2 - 3 r z (\[Phi]^2 + \[Phi] - 2) + 2 z^2 (\[Phi]^2 + \[Phi] + 1)))/(2 r^2 (\[Phi] - 1)^3)]; (* define integral function *) myIntFun[r_?NumericQ, p_?NumericQ, d_?NumericQ, \[Phi]_?NumericQ] := NIntegrate[wFun[z, r, p, d] qFun[z, r, \[Phi]], {z, 0, 10}]; (* for now, create table varying just phi from 0.01 to 0.99 and fit a \ formula to this 1D data. *) phiTable = Table[ {phi, myIntFun[1, 1, 2, phi]}, {phi, 0.01, 0.99, 0.01}]; (* find a formula for data *) theF[x_] = FindFormula[phiTable, x] (* superimpose ListPlot of points with fitted function *) lp = ListPlot[phiTable, PlotStyle -> {Black, PointSize[0.01]}] p1 = Plot[theF[x], {x, 0.01, 0.99}, PlotStyle -> Red] Show[{lp, p1}]