Error using integral function & function handles in a loop

'Outliers.m' is called from a higher level .m file. The variables ('cc', 'mastercc' etc) are all defined in the higher level file, and set as globals for access by Outliers.m. The purpose of the code is to identify outliers using Chauvenets Criterion, and for this, I have to calculate the integral of the guassian distribution, using the Integral function. The code works and gives sensible values when I enter specific variables as a test with NO loop, but I cannot get it to work inside a loop. My data set ('cc' = (7x30)) is comprised of 7 individual samples, each 1x30 , all of which need to be analyzed. I have read through the guidance on Integral and function handles, and various forums, but cannot seem to find the solution...Any help or guidance would be very much appreciated.... Here is my code:

    n = 7
    for x = 1:n
      for y = 1:30
      z(x,y) = abs((cc(x,y) - mastercc(1,y))/masterccstd(1,y));
      xmax(x,y) = mastercc(1,y)+z(x,y)*masterccstd(1,y);
      xmin(x,y) = mastercc(1,y)-z(x,y)*masterccstd(1,y);
      p(x,y) = 1/(masterccstd(1,y)*(sqrt(2*pi)));

      fun = @(z)exp(-1/2)*z.^2;
      q(x,y) = integral(fun(z(x,y)),xmin(x,y),xmax(x,y));

      pq(x,y) = p(x,y)*q(x,y); % probability
      value(x,y) = n*(1/pq(x,y));
      count(x,y) = logical(value(x,y) <0.5);
      badbins(x)=sum(count(x,:));
      end
  end

The error from this code (bp error stops at the 'integral' line) is:

"Error using integral (line 83) First input argument must be a function handle. Error in outliers (line 30) q(x,y) = integral(fun(z(x,y)),xmin(x,y),xmax(x,y));"

note: 'fun' DOES show up as a function handle in the workspace (confirmed by: isa(fun,'function_handle') = 1)

The variables at the error point are : x = 1, y = 1, z = 0.2502, xmax = 1.9428, xmin = 1.9104.

If it helps, here is the same code but with no loops, and the variables replaced with constants (to clarify - this code does work - it is the code in a loop which doesn't):

n = 7;
z = abs(4.3794 - 1)/1;
xmax = 1+z*1;
xmin = 1-z*1;
p = 1/(1*2.5066);
fun = @(z) exp(-1/2*z.^2);
q = integral(fun,xmin,xmax);
pq = p*q; value = 6*1/pq;
count = logical(value <0.5);
badbins=sum(count);

Answer 1

[Andrei Bobrov]

z = abs(bsxfun(@minus,cc,mastercc));
xmin = bsxfun(@minus,mastercc,z);
xmax = bsxfun(@plus,mastercc,z);
p = 1./(masterccstd*(sqrt(2*pi)));
fun = @(z)exp(-1/2)*z.^2;
q = arrayfun(@(ii,jj)integral(fun,ii,jj),xmin,xmax);
pq = bsxfun(@times,p,q);
value = n./pq;
count = logical(value <0.5);
badbins=sum(count,2);

[gingermonster]

@andrei, thank you very much for taking the time to answer my question. 'bsxfun' certainly helps to remove the need of loops in this code.