function [K, F] = CalculateLaplace_mult(model, lambda_wall, lambda_fluid, lambda_air)
mesh = model.Mesh;

nodes    = mesh.Nodes;
elements = mesh.Elements;
Nnodes = size(nodes,2);
Nelems = size(elements,2);

K = sparse(Nnodes, Nnodes);
F = zeros(Nnodes,1);

% --- build region array manually ---
regions = zeros(Nelems,1);
wallElems  = findElements(mesh,'region','Face',1);
fluidElems = findElements(mesh,'region','Face',2);
airElems   = findElements(mesh,'region','Face',3);
regions(wallElems)  = 1;
regions(fluidElems) = 2;
regions(airElems)   = 3;

for e = 1:Nelems
    vert = elements(:,e);
    x = nodes(1,vert);
    y = nodes(2,vert);

    Ae = polyarea(x,y);

    b = [y(2)-y(3); y(3)-y(1); y(1)-y(2)];
    c = [x(3)-x(2); x(1)-x(3); x(2)-x(1)];

    % conductivity by material
    switch regions(e)
        case 1
            lambda = lambda_wall;
        case 2
            lambda = lambda_fluid;
        case 3
            lambda = lambda_air;
        otherwise
            error('Element %d has no material assignment', e);
    end
    Ke = (lambda/(4*Ae)) * (b*b.' + c*c.');
    K(vert,vert) = K(vert,vert) + Ke;
end
end