diff --git a/Project/CalculateLaplace_mult.m b/Project/CalculateLaplace_mult.m
new file mode 100644
index 0000000..ff274ba
--- /dev/null
+++ b/Project/CalculateLaplace_mult.m
@@ -0,0 +1,45 @@
+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