diff --git a/Project/Main.m b/Project/Main.m
new file mode 100644
index 0000000..fa32fe4
--- /dev/null
+++ b/Project/Main.m
@@ -0,0 +1,281 @@
+% Axisymmetric mug 
+% r–z plane, rotation around r = 0
+clear; clc; close all;
+
+%% General values that we use in the entire script 
+%Task 3) Thermal conductivity 
+lambda_wall  = 30;    % ceramic 
+lambda_fluid = 0.6089;    % water 
+lambda_air   = 0.026;  % air
+%Task 4) 
+alpha = 10;      % heat transfer coefficient
+u_out = 18;      % ambient air temperature (°C)
+%Task 6) Volumetric heat capacities: rho * c_p
+% Densities [kg/m^3]
+rho_wall  = 2400;   % ceramic
+rho_fluid = 1000;   % water
+rho_air   = 1.2;    % air
+% Specific heats [J/(kg*K)]
+cp_wall  = 900;
+cp_fluid = 4184;
+cp_air   = 1005;
+% Volumetric heat capacities [J/(m^3*K)]
+c_wall  = rho_wall  * cp_wall;
+c_fluid = rho_fluid * cp_fluid;
+c_air   = rho_air   * cp_air;
+
+%% Task 1: Mesh definition
+%  Create PDE model
+model = createpde();
+
+%  Points (meters)
+A = [0,     0];         B = [0.055, 0];
+C = [0.083, 0.105];     H = [0.078, 0.105];
+F = [0.050, 0.005];     E = [0,     0.005];
+G = [0.067, 0.066];     I = [0,     0.066];
+D = [0,     0.105];
+
+%  Geometry matrix (edges)
+g1  = [2; A(1); E(1); A(2); E(2); 1; 0]; % Axis - ceramic
+g2  = [2; E(1); I(1); E(2); I(2); 2; 0]; % Axis - fluid
+g3  = [2; I(1); D(1); I(2); D(2); 3; 0]; % Axis - air 
+g4  = [2; A(1); B(1); A(2); B(2); 1; 0]; % Outer ceramic
+g5  = [2; B(1); C(1); B(2); C(2); 1; 0]; % Outer ceramic
+g6  = [2; C(1); H(1); C(2); H(2); 1; 3]; % Top rim: (C-H) ceramic-air
+g7  = [2; H(1); F(1); H(2); F(2); 1; 3]; % Inner ceramic wall (H-F) ceramic-air
+g8  = [2; F(1); E(1); F(2); E(2); 1; 2]; % Inner ceramic bottom (F-E) ceramic-fluid
+g9  = [2; F(1); G(1); F(2); G(2); 2; 3]; % Fluid surface (F-G) fluid-air
+g10 = [2; G(1); I(1); G(2); I(2); 2; 3]; % Fluid surface (G-I) fluid-air
+g11 = [2; D(1); H(1); D(2); H(2); 3; 0]; % Air top boundary: (D-H) air-outside
+
+% Assemble geometry
+g = [g1 g2 g3 g4 g5 g6 g7 g8 g9 g10 g11];
+geometryFromEdges(model, g);
+
+% figure(1);
+% pdegplot(model, 'EdgeLabels','on', 'FaceLabels','on');
+% axis equal;
+% title('Geometry with edge and face labels');
+
+%  Generate mesh, linear and 3 nodes per element
+mesh = generateMesh(model, 'Hmax', 0.002, 'GeometricOrder','linear');
+
+% figure(2);
+% pdemesh(model);
+% axis equal;
+% title('Generated mesh');
+
+%% Task 2: Direct solver with constant lambda 
+nodes = mesh.Nodes;         % coordinates of all mesh nodes
+elements = mesh.Elements;   % which nodes make up each triangle element
+
+Nnodes = size(nodes,2);
+Nelems = size(elements,2);
+
+% Define material properties (for simplicity, lambda = 1 everywhere)
+lambda = ones(Nelems,1);    % thermal conductivity
+
+% Initialize global stiffness matrix and RHS
+K = sparse(Nnodes, Nnodes);
+F = zeros(Nnodes,1);
+
+% Assemble K and F
+for e = 1:Nelems %Loop over each triangle element 
+    vert = elements(:,e);   %nodes of element
+    x = nodes(1,vert);
+    y = nodes(2,vert);
+
+    Ae = polyarea(x,y);     % Compute area of the triangle 
+
+    % Linear triangle gradients
+    b = [y(2)-y(3); y(3)-y(1); y(1)-y(2)]; % derivative with respect to x 
+    c = [x(3)-x(2); x(1)-x(3); x(2)-x(1)]; % derivative with respect to y
+
+    % Element stiffness matrix
+    Ke = (lambda(e)/(4*Ae)) * (b*b.' + c*c.');
+
+    % Assemble
+    K(vert,vert) = K(vert,vert) + Ke;
+
+    % Element load vector (f=0)
+    F(vert) = F(vert) + zeros(3,1);
+end
+
+% Find boundary nodes (couldn't find a better way)
+edgesAll = [elements([1 2],:), elements([2 3],:), elements([3 1],:)]; % all edges -> 2x(3*Nelems) awway since each triangle has 3 edges 
+edgesSorted = sort(edgesAll,1); % sort nodes of each edge, ensure that [i,j] and [j,i] are the same
+[~,~,ic] = unique(edgesSorted','rows'); % identifies unique edges and assignes indices ic 
+counts = accumarray(ic,1); % counts how many times each edge appears in the mesh
+boundaryEdges = find(counts==1); %edges belonging to only 1 element: hence on the boundary 
+boundaryNodes = unique(edgesSorted(:,boundaryEdges)); % nodes belonging to these boundary edges 
+
+% Direct solver
+% Enforce Dirichlet BC strongly
+K(boundaryNodes,:) = 0;
+K(:,boundaryNodes) = 0;
+K(boundaryNodes,boundaryNodes) = speye(length(boundaryNodes));
+F(boundaryNodes) = 0;
+
+% Direct solve
+u = K \ F;
+
+% figure(3)
+% pdeplot(model, 'XYData', u, 'Mesh','on');
+% axis equal;
+% title('Stationary Dirichlet solution');
+% colorbar;
+
+%% Task 3: Laplace with multiple lambdas
+[K, F] = CalculateLaplace_mult(model, lambda_wall, lambda_fluid, lambda_air);
+
+% Enforce Dirichlet BC strongly
+K(boundaryNodes,:) = 0;
+K(:,boundaryNodes) = 0;
+K(boundaryNodes,boundaryNodes) = speye(length(boundaryNodes));
+F(boundaryNodes) = 0;
+
+% Direct solve
+u = K \ F;
+
+% % Plot solution
+% figure(4)
+% pdeplot(model,'XYData',u, 'Mesh','on');
+% axis equal
+% title('Task 3: Stationary solution with multiple conductivities');
+% colorbar
+
+%% Task 4: Robin boundary condition 
+[K, F] = CalculateLaplace_mult(model, lambda_wall, lambda_fluid, lambda_air);
+[K, F] = ApplyRobinBC_mult(model, K, F, alpha, u_out);
+
+% Direct solve
+u = K \ F;
+
+% figure(5)
+% pdeplot(model, 'XYData', u, 'Mesh','on');
+% axis equal
+% title('Task 4: Stationary solution with Robin BC');
+% colorbar
+
+%% Task 5: Axisymmetric Laplace + Robin BC 
+[K, F] = CalculateLaplace_mult_rot(model, lambda_wall, lambda_fluid, lambda_air);
+[K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out);
+
+% Direct solve
+u = K \ F;
+
+% figure(6)
+% pdeplot(model, 'XYData', u, 'Mesh','on');
+% axis equal
+% title('Task 5: Axisymmetric stationary solution with Robin BC');
+% colorbar
+
+%To see it in 3D paste here the code in "AdditionalPlotCodes.txt". 
+
+%% Task 6: Axisymmetric mass matrix
+M = sparse(Nnodes, Nnodes);
+M = AddMass_mult_rot(model, M, c_wall, c_fluid, c_air);
+
+%% Task 7: Initial solution
+u0 = Init_Solution_mult(model, 18, 80, 18);
+
+% figure(7)
+% pdeplot(model, 'XYData', u0, 'Mesh','on');
+% axis equal
+% title('Initial temperature distribution');
+% colorbar
+
+%% Task 8: Time-dependent simulation (explicit scheme)
+tau = 0.5;              % time step in seconds 
+T_end = 400;            % total simulation time (seconds)
+Nt = ceil(T_end/tau);   % number of time steps
+
+A = (1/tau)*M+K; % Left-hand side matrix
+u = u0;
+
+for k = 1:Nt
+    b = (1/tau)*M*u + F;   % F is the load vector, F=0 
+    u_next = A\b;
+    u = u_next; % Update
+    
+    if mod(k,20) == 0
+        % figure(8)
+        % pdeplot(model, 'XYData', u, 'Mesh','on');
+        % axis equal
+        % title(['Temperature at t = ', num2str(k*tau), ' s']);
+        % colorbar
+        % drawnow
+    end
+end
+
+%To see the 9 snapshots paste here the codes in "AdditionalPlotCodes.txt"
+
+%% Task 9 (i): Heating time using inner ceramic wall temperature
+T_target = 67;   % [°C]
+[K, F] = CalculateLaplace_mult_rot(model, lambda_wall, lambda_fluid, lambda_air);
+[K, F] = ApplyRobinBC_mult_rot(model, K, F, alpha, u_out);
+A = (1/tau)*M+K; % Left-hand side matrix
+
+innerWallNodes = findNodes(model.Mesh,'region','Edge',8); % Edge 8 = ceramic–fluid interface 
+u = u0;
+
+% Storage
+timeVec = (0:Nt-1)' * tau;
+innerWallTemp = zeros(Nt,1);
+Twarm = NaN;
+
+for k = 1:Nt
+    b = (1/tau)*M*u + F;
+    
+    u = A\b;
+
+    % Average inner wall temperature
+    innerWallTemp(k) = mean(u(innerWallNodes));
+
+    % Check heating criterion
+    if innerWallTemp(k) >= T_target
+        Twarm = k * tau;
+        fprintf('Task 9 (i): Inner wall reaches %.1f°C at T = %.1f s\n', ...
+                T_target, Twarm);
+        break
+    end
+end
+
+% figure(9)
+% plot(timeVec(1:k), innerWallTemp(1:k), 'LineWidth', 2)
+% hold on
+% yline(T_target,'r--','67°C','LineWidth',1.5)
+% xlabel('Time [s]')
+% ylabel('Average inner wall temperature [°C]')
+% title('Heating of the inner ceramic wall')
+% grid on
+
+%% CHECK: Insulated mug – transient redistribution
+[K_neu, F_neu] = CalculateLaplace_mult_rot(model,lambda_wall,lambda_fluid,lambda_air);
+
+M_neu = sparse(Nnodes, Nnodes);
+M_neu = AddMass_mult_rot(model,M_neu,c_wall,c_fluid,c_air);
+
+% Time stepping parameters
+tau = 0.5;
+T_end = 400;
+Nt = ceil(T_end / tau);
+
+A_neu = (1/tau) * M_neu + K_neu;
+
+% Initial condition
+u = u0;
+
+for k = 1:Nt
+    b = (1/tau) * M_neu * u + F_neu;
+    u = A_neu \ b;
+
+    if mod(k,20) == 0
+        % figure(10)
+        % pdeplot(model, 'XYData', u, 'Mesh','on');
+        % axis equal
+        % title(['Insulated mug, t = ', num2str(k*tau), ' s']);
+        % colorbar
+        % drawnow
+    end
+end
\ No newline at end of file