How do I write this double integral function in MATLAB and create it as a separate function file?

This is the double integral Double integral snapshot

This is so far what I've done Double integral attempt

Seems like such a simple thing. I know. I have bad teachers. Thanks for the help! P.S the version I'm using R2013a

**Edit: Let me just copy and paste the code

function [] = rho_s_fun()
%UNTITLED4 Summary of this function goes here
% Amistupid.docx
%   Detailed explanation goes here
    fun =@(x,y) x*y*(x^2+y^2+25)^3/2;
    Q = integral2(fun,0,0,1,1);

end

the error is:

>> rho_s_fun
Error using  * 
Inner matrix dimensions must agree.

Error in rho_s_fun>@(x,y)x*y*(x^2+y^2+25)^3/2 (line
5)
    fun =@(x,y) x*y*(x^2+y^2+25)^3/2;

Error in integral2Calc>integral2t/tensor (line 238)
            Z1 = FUN(X(VTSTIDX),Y(VTSTIDX)); NFE =
            NFE + 1;

Error in integral2Calc>integral2t (line 56)
[Qsub,esub] = tensor(thetaL,thetaR,phiB,phiT);

Error in integral2Calc (line 10)
    [q,errbnd] =
    integral2t(fun,xmin,xmax,ymin,ymax,optionstruct);
    
Error in integral2 (line 107)
    Q =
    integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);
    
Error in rho_s_fun (line 6)
    Q = integral2(fun,0,0,1,1); 

The main code where I'm trying to link the double integral(not sure if link is the proper word) is:

clear
format short

%The following part is for part (a)
dx = input('Enter the integration increment dx... \n >');
dy = input('Enter the integration increment dy... \n >');
total_charge=0;
for x=0:dx:1, 
for y=0:dy:1,
total_charge=rho_s_fun(x,y)*dx*dy+total_charge;
end
end
disp(sprintf('The total charge computed by manual double'));
disp(sprintf(' integration is %d nC',total_charge));
total_charge=dblquad(@rho_s_fun, 0, 1, 0, 1);
disp(sprintf('The total charge computed by the functional'))
disp(sprintf(' double integration is %d n',total_charge))

%The following for part (b)
Etotal=0;
for x=0:dx:1,
for y=0:dx:1,
rminusrprime=[-x -y 5];
numerator=rho_s_fun(x,y)*1e-9*dx*dy*rminusrprime;
denominator=4*pi*1e-9/(36*pi)*norm(rminusrprime)^3;
Etotal=Etotral+numerator/denominator;
end
end
disp(sprintf('The electric field computed by manual double'));
disp(sprintf(' integration is (%d, %d, %d) V/m',... 
    Etotal(1), Etotal(2), Etotal(1)))