Eshelby's tensor for three dimensional mesh

Dear all,

Seasons' greetings.

I am trying to find out the Eshelby's tensor for superellipsoid system. I am following the protocol of finding flux integral from this link https://www3.nd.edu/~nancy/Math20550/Matlab/Assignments/SurfaceIntegrals/surfaceintegrals.html#27

The code returns me the symbolic coefficients, rather than any numeric value.

Please help.

%%%%%%%%%%%%%%% designing superellipsoid %%%%%%%%%%%%%%%%%
                 %declare constants%
a1=25; a2=25; a3=25; epsilon1=1; epsilon2=1; epsilon3=1;
n=100;
  
                 %declare variables%
 etamax=pi/2; etamin=-pi/2;
 wmax=pi; wmin=-pi;
 deta=(etamax-etamin)/n;
 dw=(wmax-wmin)/n;
 [i,j] = meshgrid(1:n+1,1:n+1);
 eta = etamin + (i-1) * deta; w   = wmin + (j-1) * dw;
  
%%%%%%%%%%%%%%% crating symbolic variables %%%%%%%%%%%%%%%%%
syms eta w

ellip=[(a1.*sign(cos(eta)).*abs(cos(eta)).^epsilon1.*sign(cos(w)).*abs(cos(w)).^epsilon1),(a2.*sign(cos(eta)).*abs(cos(eta)).^epsilon2.*sign(sin(w)).*abs(sin(w)).^epsilon2),(a3.*sign(sin(eta)).*abs(sin(eta)).^epsilon3)];
F=[(a1.*sin(eta).*cos(w)) (a2.*sin(eta).*sin(w)) (a3.*cos(eta))];
nds=simplify(cross(diff(ellip,eta),diff(ellip,w)));
Fpar=subs(F,[(a1.*sin(eta).*cos(w)) (a2.*sin(eta).*sin(w)) (a3.*cos(eta))],ellip);
Fj=@(Fpar,nds)Fpar.*transpose(nds);
flux=(symint2(Fj,eta,-pi,pi,w,-(pi/2),(pi/2)));