Calculating the posterior distribution in multivariate Gaussian processes
Assuming model.KuyDinv{k,q} and model.Kyu{q,r} are a matrices (as the name would imply), there is little you can do. You can move the initialization of KuyDinvKyu outside the loop, and eliminate the branch to compute the diagonals: 

% define it once -- saves many calls to 'zeros'
KuyDinvKyu_0 = zeros(model.k);

for k = 1:model.nlf
    for r = 1:model.nlf

        KuyDinvKyu = KuyDinvKyu_0;

        for q = 1:model.nout
            KuyDinvKyu = KuyDinvKyu + model.KuyDinv{k,q}*model.Kyu{q,r};
        end

        % removed IF      
        model.A{k,r} = KuyDinvKyu;

    end
end

% without IF
for c = 1:model.nlf    
    model.A{c,c} = model.A{c,c} + model.Kuu{c};    
end


If model.KuyDinv{k,q} and/or model.Kyu{q,r} contain scalars, well then we can optimize this much further for sure. So I need to know the size and type of data that model.Kyu{q,r} and  model.KuyDinv{k,q} will contain. 

It might also be that your overall data design (e.g., the choice to use cells) is flawed and causes inefficiencies. So it could also be helpful to see some more surrounding code, so I can determine if there is some improvement to be made in that respect.