Slow performance using closures

I'm coding a solution for Poisson equation on a 2d rectangle using finite elements. In order to simplify the code I store handles to the basis functions in an array and then loop over these basis functions to create my matrix and right hand side. The problem with this is that even for very coarse grids it is prohibitively slow. For a 9x9 grid (using Dirichlet BC, there are 49 nodes to solve for) it takes around 20 seconds. Using the profile I've noticed that around half the time is spent accessing (not executing) my basis functions.

The profiler says:

 matrix_assembly>@(x,y)bilinearBasisFunction(x,y,xc(k-1),xc(k),xc(k+1),yc(j-1),yc(j),yc(j+1)) (156800 calls, 11.558 sec)

The self time (not executing the bilinear basis code) is over 9 seconds. Any ideas as to why this might be so slow?

Here's some of the code (matrix_assembly.m), I can post more if needed:

%%setting up the basis functions, storing them in cell array
basisFunctions = cell(nu, 1); %nu is #unknowns 
i = 1;
for j = 2:length(yc) - 1
for k = 2:length(xc) - 1
  basisFunctions{i} = @(x,y) bilinearBasisFunction(x,y, xc(k-1), xc(k),...
      xc(k+1), yc(j-1), yc(j), yc(j+1)); %my code for bilinear basis functions
  i = i+1;
end
end

 %%Assemble matrices and RHS
M = zeros(nu,nu);
S = zeros(nu,nu);
F = zeros(nu, 1);

 for iE = 1:ne %ne is # elements (64 is 9x9 case)
 for iBF = 1:nu %nu is # unknowns (49 in 9x9 case)
    [z1, dx1, dy1] = basisFunctions{iBF}(qx(iE), qy(iE));

    F(iBF) = F(iBF) + z1*forcing_handle(qx(iE),qy(iE))/ae(iE);

    for jBF = 1:nu
        [z2, dx2, dy2] = basisFunctions{jBF}(qx(iE), qy(iE));

        %M(iBF,jBF) = M(iBF,jBF) + z1*z2/ae(iE);
        S(iBF,jBF) = S(iBF, jBF) + (dx1*dx2 + dy1*dy2)/ae(iE);
    end        
 end
 end