Бавна производителност на Matlab при използване на затваряния

Кодирам решение за уравнение на Поасон върху 2d правоъгълник, използвайки крайни елементи. За да опростя кода, съхранявам манипулатори на базовите функции в масив и след това преминавам през тези базови функции, за да създам моята матрица и дясната страна. Проблемът с това е, че дори за много груби мрежи е непосилно бавен. За мрежа 9x9 (използвайки Dirichlet BC, има 49 възела за решаване) отнема около 20 секунди. Използвайки профила, забелязах, че около половината от времето се изразходва за достъп (не изпълнение) на моите основни функции.

Профайлърът казва 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), самото време (без изпълнение на кода на билинейната основа) е над 9 секунди. Някакви идеи защо това може да е толкова бавно?

Ето част от кода, мога да публикувам още, ако е необходимо:

%% 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
for iBF = 1:nu
    [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