我正在尝试解决氦原子量子力学中的物理问题，该解决方案需要对 2 个变量进行数值积分。但是，当我尝试运行下一个代码时

e = 1.6e-19;
eps0 = 8.85e-12;
a0 = 0.53e-10;
c = e^2/(4*pi*eps0);
psi1 =@(r) ((2^1.5)/(sqrt(pi)*a0^1.5))*exp(-2*r/a0);
psi2 =@(r) ((2^1.5)/(4*sqrt(2*pi)*a0^1.5))*(2-r/a0)*exp(-r/a0);
I =@(r1,r2) c*(psi2(r1).*psi1(r2).*(1./abs(r1-r2)).*psi2(r1).*psi1(r2));
J =@(r1,r2) c*(psi2(r1).*psi1(r2).*(1./abs(r1-r2)).*psi1(r1).*psi2et(r2));
integral2(I,0,inf,0,inf)
integral2(J,0,inf,0,inf)

matlab 说

Error using  * 
Inner matrix dimensions must agree.

Error in qmex7>@(r)((2^1.5)/(4*sqrt(2*pi)*a0^1.5))*(2-r/a0)*exp(-r/a0)

Error in qmex7>@(r1,r2)c*(psi2(r1).*psi1(r2).*(1./abs(r1-r2)).*psi2(r1).*psi1(r2))

Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)
    @(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), ...

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

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

Error in integralCalc (line 83)
        [q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions)

Error in integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x)) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integralCalc/iterateScalarValued (line 314)
                fx = FUN(t);

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

Error in integralCalc (line 83)
        [q,errbnd] = vadapt(@AToInfInvTransform,interval);

Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

Error in integral2Calc (line 7)
    [q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 106)
    Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

知道为什么吗？