四重积分的相关问题

[zqggsimida]

各位好，小弟因为学习需要刚接触matlab，遇到了一个问题，因为身边没有学习matlab的人，特此前来求助。就是我需要计算图片中的四重积分，因为matlab好像没有直接求解四重积分的函数，所以我采用了一重符号加三重数值的方法如下面的代码所示(版本为2016)
syms x1 x2 y1 y2
L=sqrt((x1-x2).^2+(y1-y2).^2+1);
I1=int(1/(pi*L.^4),y2,-sqrt(1-x2.^2),sqrt(1-x2.^2));
fun1=eval(['@(x1,y1,x2)',vectorize(I1)]);
x1min=-1;x1max=1;
y1min=@(x1)-sqrt(1-x1.^2);y1max=@(x1)sqrt(1-x1.^2);
x2min=-1;x2max=1;
y=integral3(fun1,x1min,x1max,y1min,y1max,x2min,x2max)

结果运行的时候出了如下错误
出错 integral2Calc>integral2t/tensor (line 228)
Z = FUN(X,Y); NFE = NFE + 1;

出错 integral2Calc>integral2t (line 55)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

出错 integral2Calc (line 9)
[q,errbnd] = integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);

出错 integral3/innerintegral (line 137)
Q1 = integral2Calc( ...

出错 integralCalc/iterateScalarValued (line 314)
fx = FUN(t);

出错 integralCalc/vadapt (line 132)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen);

出错 integralCalc (line 75)
[q,errbnd] = vadapt(@AtoBInvTransform,interval);

出错 integral3 (line 121)
Q = integralCalc(@innerintegral,xmin,xmax,integralOptions);

所以想请教一下各位我的代码问题所在（因为刚学,可能写的菜，还请各位担待），不甚感激！！感谢各位！！