Numerical Integration in MatLab -

i trying compute integral using matlab's dblquad. first wrote script function. code is

function z = intrect(x, y) %the variables z0=20; x0=15; y0=20; rl=sqrt(x.^2+y.^2+(z0/2)^2); theta=acos(z0./(2*rl)); phi=atan(y./x); %the function z=(x0-z0*tan(theta).*cos(phi))*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; z=z/rl.^3; 

to compute numerical integral type in command window

z0=20;x0=15;y0=20; q = dblquad(@intrect,0,x0/2,0,y0/2,1.e-6); 

i error saying that

??? error using ==> mtimes inner matrix dimensions must agree.  error in ==> intrect @ 8 z=(x0-z0*tan(theta).*cos(phi))*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4;  error in ==> quad @ 77 y = f(x, varargin{:});  error in ==> dblquad>innerintegral @ 84     q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), varargin{:});  error in ==> quad @ 77 y = f(x, varargin{:});  error in ==> dblquad @ 60 q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ... 

what doing wrong that?

edit

replacing

z=(x0-z0*tan(theta).*cos(phi))*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; 

with

z=(x0-z0*tan(theta).*cos(phi)).*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; 

i new error

??? index exceeds matrix dimensions.  error in ==> quad @ 85 if ~isfinite(y(7))  error in ==> dblquad>innerintegral @ 84     q(i) = quadf(intfcn, xmin, xmax, tol, trace, y(i), varargin{:});  error in ==> quad @ 77 y = f(x, varargin{:});  error in ==> dblquad @ 60 q = quadf(@innerintegral, ymin, ymax, tol, trace, intfcn, ... 

as dblquad says, input x vector , y scalar , output z vector. in integrand function function of x vector (e.g., rl) , you'll need careful use element-wise operators appropriate. did not last line.

also, consider passing initial value parameters via function handle rather duplicating them inside integrand function:

function dblquadtest z0 = 20; x0 = 15; y0 = 20; f = @(x,y)intrect(x,y,x0,y0,z0); q = dblquad(f,0,x0/2,0,y0/2); % 1e-6 default tolerance  function z = intrect(x, y, x0, y0, z0) %the variables rl=sqrt(x.^2+y.^2+(z0/2)^2); theta=acos(z0./(2*rl)); phi=atan(y./x); %the function z=(x0-z0*tan(theta).*cos(phi)).*(y0-z0*tan(theta).*sin(phi))*(z0/2)^4; z=z./rl.^3;