diff --git a/Driven_Cavity_Solver.m b/Driven_Cavity_Solver.m
new file mode 100644
index 0000000000000000000000000000000000000000..6c17e5e9d0e6ef83517a07eff2cbe3c7595dda13
--- /dev/null
+++ b/Driven_Cavity_Solver.m
@@ -0,0 +1,283 @@
+function Driven_Cavity_Solver()
+ %% Driven Cavity, Finite Volume Method
+ % Paul Bataillie - paul.bataillie@student.isae-supaero.fr
+ % ---------------------------------------------------------------------
+
+ close all
+ clc
+ clear
+
+ %% Constants and parameters
+ Lx = 1;
+ Ly = 1;
+ rho = 1000;
+ mu = 1e-3;
+
+ % Underrelaxation factors
+ e = 0.5;
+ ep = 0.2;
+
+ % Max iterations
+ numItr = 10000;
+ numItrP = 800;
+ dispResiduals = 20;
+
+ % Number of cells
+ Npx = 20;
+
+ disp(append("INFO: Initializing Case ",num2str(Npx),"x",num2str(Npx),"."))
+ disp(append("INFO: Under-Relaxation Factors | e = ",num2str(e)," | ep = ",num2str(ep)))
+
+ Npy = Npx; % The pressure matrix (grid) is a square
+
+ Nvx = Npx+2; % The v-matrix contains two more columns and one more row
+ Nvy = Npx+1; % than the pressure matrix (including dummy cells)
+
+ Nux = Npx+1; % The u-matrix contains one more column and two more rows
+ Nuy = Npy+2; % than the pressure matrix (including dummy cells)
+
+ % Cell size
+ dx = Lx/(Npx);
+ dy = Ly/(Npy);
+
+ % Residual tolerances
+ RCE = 1e-5;
+ RME = 1e-7;
+
+ %Reynolds number
+ Re = 1000;
+ Re1 = 1/Re;
+
+ % Some variables
+ Dew = Re1*dy/dx;
+ Dns = Re1*dx/dy;
+
+ % Initialization
+ u = zeros(Nuy,Nux);
+ unew = zeros(Nuy,Nux);
+
+ v = zeros(Nvy,Nvx);
+ vnew = zeros(Nvy,Nvx);
+
+ p = zeros(Npy,Npx);
+ pp = zeros(Npy+2,Npx+2);
+ pp1 = zeros(Npy+2,Npx+2);
+
+ auP = 0*ones(Nuy,Nux);
+ avP = 0*ones(Nvy,Nvx);
+
+ uo = u;
+ vo = v;
+
+ scrsz = get(0,'ScreenSize');
+
+ %Boundray Conditions
+ unew(end,:) = 2 - unew(end-1,:);
+ u(end,:) = 2 - u(end-1,:);
+
+ %Residuals Variables
+ continuity_res_scale_value = 0;
+ u_residual_vector = zeros(numItr,1);
+ v_residual_vector = zeros(numItr,1);
+ c_residual_vector = zeros(numItr,1);
+
+ %% Start of the solution process
+ format long
+
+ disp("INFO: Solving...")
+
+ for i=1:numItr
+ %% Solve momentum equations
+ %% ------
+ % Mex
+ % Average velocities for the flux terms
+ uue = 0.5*(u(2:end-1,3:end) + u(2:end-1,2:end-1));
+ uuw = 0.5*(u(2:end-1,1:end-2) + u(2:end-1,2:end-1));
+ vun = 0.5*(v(2:end,2:end-2) + v(2:end,3:end-1));
+ vus = 0.5*(v(1:end-1,2:end-2) + v(1:end-1,3:end-1));
+
+ %Flux terms
+ Fue = dy*uue;
+ Fuw = dy*uuw;
+ Fun = dx*vun;
+ Fus = dx*vus;
+
+ % Coefficients
+ auE = Dew + max(0,-Fue);
+ auW = Dew + max(0,Fuw);
+ auN = Dns + max(0,-Fun);
+ auS = Dns + max(0,Fus);
+ auP = auE + auW + auN + auS + Fue + Fun - Fuw - Fus;
+
+ % Solve MEx
+ unew(2:end-1,2:end-1) = (auE.*u(2:end-1,3:end) + auW.*u(2:end-1,1:end-2) + auN.*u(3:end,2:end-1) + auS.*u(1:end-2,2:end-1) - (p(1:end,2:end) - p(1:end,1:end-1))*dy) ./ (auP);
+
+ % Update boundary conditions
+ unew(end,:) = 2 - unew(end-1,:);
+ unew(:,1) = 0;
+ unew(:,end) = 0;
+ unew(1,:) = - unew(2,:);
+
+ %% ------
+ % Mey
+ % Average velocities for the flux terms
+ uve = 0.5*(u(2:end-2,2:end) + u(3:end-1,2:end));
+ uvw = 0.5*(u(2:end-2,1:end-1) + u(3:end-1,1:end-1));
+ vvn = 0.5*(v(2:end-1,2:end-1) + v(3:end,2:end-1));
+ vvs = 0.5*(v(1:end-2,2:end-1) + v(2:end-1,2:end-1));
+
+ %Flux terms
+ Fve = dy*uve;
+ Fvw = dy*uvw;
+ Fvn = dx*vvn;
+ Fvs = dx*vvs;
+
+ % Coefficients
+ avE = Dew + max(0,-Fve);
+ avW = Dew + max(0,Fvw);
+ avN = Dns + max(0,-Fvn);
+ avS = Dns + max(0,Fvs);
+ avP = avE + avW + avN + avS + Fve + Fvn - Fvw - Fvs;
+
+ % Solve MEy
+ vnew(2:end-1,2:end-1) = (avE.*v(2:end-1,3:end) + avW.*v(2:end-1,1:end-2) + avN.*v(3:end,2:end-1) + avS.*v(1:end-2,2:end-1) - (p(2:end,1:end) - p(1:end-1,1:end))*dx) ./ (avP);
+
+ % Update boundary conditions
+ vnew(1,:) = 0;
+ vnew(end,:) = 0;
+ vnew(:,1) = - vnew(:,2);
+ vnew(:,end) = - vnew(:,end-1);
+
+ %% ------
+ % Pressure correction
+ pp1 = zeros(Npy+2,Npx+2);
+
+ apE = zeros(Npy,Npx);
+ apN = zeros(Npy,Npx);
+ apW = zeros(Npy,Npx);
+ apS = zeros(Npy,Npx);
+
+ apE(:,1:end-1) = dy^2./auP;
+ apE(:,end) = 0;
+
+ apN(1:end-1,:) = dx^2./avP;
+ apN(end,:) = 0;
+
+ apW(:,2:end) = dy^2./auP;
+ apW(:,1) = 0;
+
+ apS(2:end,:) = dx^2./avP;
+ apS(1,:) = 0;
+
+ ap = apE + apW + apN + apS;
+ bp = (-unew(2:end-1,2:end)+unew(2:end-1,1:end-1))*dy + (-vnew(2:end,2:end-1)+vnew(1:end-1,2:end-1))*dx;
+
+ % Loop for p'
+ for j=1:numItrP
+ pp1(2:end-1,2:end-1) = ( apE.*pp(2:end-1,3:end) + apW.*pp(2:end-1,1:end-2) + apN.*pp(3:end,2:end-1) + apS.*pp(1:end-2,2:end-1) + bp ) ./ ap;
+
+ % Underrelaxation
+ pp = ep*pp1 + (1-ep)*pp;
+
+ end %...of p' loop
+
+ % Update p
+ p = p + pp(2:end-1,2:end-1);
+
+ % Update u and v in internal points
+ u(2:end-1,2:end-1) = e * ( unew(2:end-1,2:end-1) + ( (dy) * (-pp(2:end-1,3:end-1) + pp(2:end-1,2:end-2)) )./(auP) ) + (1-e)*u(2:end-1,2:end-1);
+ v(2:end-1,2:end-1) = e * ( vnew(2:end-1,2:end-1) + ( (dx) * (-pp(3:end-1,2:end-1) + pp(2:end-2,2:end-1)) )./(avP) ) + (1-e)*v(2:end-1,2:end-1);
+
+ % Update boundary conditions
+ u(end,:) = 2 - unew(end-1,:);
+ u(:,1) = 0;
+ u(:,end) = 0;
+ u(1,:) = - unew(2,:);
+
+ v(1,:) = 0;
+ v(end,:) = 0;
+ v(:,1) = - vnew(:,2);
+ v(:,end) = - vnew(:,end-1);
+
+ if (mod(i,dispResiduals)==0)||(i < 6)
+ %% Residuals
+ % U Momentum Sclaed Residual
+ num = abs(unew(2:end-1,2:end-1).*(auP) - ( (auE.*u(2:end-1,3:end) + auW.*u(2:end-1,1:end-2) + auN.*u(3:end,2:end-1) + auS.*u(1:end-2,2:end-1) - (p(1:end,2:end) - p(1:end,1:end-1))*dy) ));
+ den = abs(unew(2:end-1,2:end-1).*(auP));
+ r_u = sum(num,"all")/sum(den,"all");
+
+ % V Momentum Scaled Residual
+ num = abs(vnew(2:end-1,2:end-1).*(avP) - ( (avE.*v(2:end-1,3:end) + avW.*v(2:end-1,1:end-2) + avN.*v(3:end,2:end-1) + avS.*v(1:end-2,2:end-1) - (p(2:end,1:end) - p(1:end-1,1:end))*dx) ));
+ den = abs(vnew(2:end-1,2:end-1).*(avP));
+ r_v = sum(num,"all")/sum(den,"all");
+
+ % Continuity Scaled Residual (scaled as in Fluent, the scaled value is the largest absolute value of the continuity residual in the first five iterations)
+ num = sum( abs( (u(2:end-1,2:end)-u(2:end-1,1:end-1))*dy + (v(2:end,2:end-1)-v(1:end-1,2:end-1))*dx ), "all");
+ if i < 6
+ continuity_res_scale_value = max(continuity_res_scale_value,num);
+ end
+ den = continuity_res_scale_value;
+ r_c = num/den;
+
+ % Monitor points
+ disp('----')
+ disp(append("numItr = ",num2str(i)," | U Residual = ",num2str(r_u)," | V Residual = ",num2str(r_v)," | C Residual = ",num2str(r_c)))
+ disp(append("Monitor Point | u = ",num2str(u(floor(Nux/4),floor(Nuy/4)))," | v = ",num2str(v(floor(Nux/4),floor(Nuy/4)))))
+
+ u_residual_vector(i)=r_u;
+ v_residual_vector(i)=r_v;
+ c_residual_vector(i)=r_c;
+
+ if (r_u<RME) && (r_v<RME) && (r_c<RCE)
+ break
+ end
+
+ end
+
+ end
+
+ %% Plotting
+
+ % Useful commands: meshgrid, streamslice, quiver, contour
+
+ u_plot = zeros(Npy+2,Npx+2);
+ v_plot = zeros(Npy+2,Npx+2);
+ u_plot(2:end-1,2:end-1) = 0.5*(u(2:end-1,1:end-1) + u(2:end-1,2:end));
+ v_plot(2:end-1,2:end-1) = 0.5*(v(1:end-1,2:end-1) + v(2:end,2:end-1));
+
+ u_plot(1,:) = 0;
+ u_plot(end,:) = 1;
+ u_plot(:,1) = 0;
+ u_plot(:,end) = 0;
+
+ v_plot(1,:) = 0;
+ v_plot(end,:) = 0;
+ v_plot(:,1) = 0;
+ v_plot(:,end) = 0;
+
+ x = zeros(1,Npx+2);
+ y = zeros(1,Npy+2);
+
+ indx = (0:1/Npx:1) + 1/(Npx*2);
+ indy = (0:1/Npy:1) + 1/(Npy*2);
+
+ x(1,2:end-1) = indx(1,1:end-1);
+ y(1,2:end-1) = indy(1,1:end-1);
+ x(1,end) = 1;
+ y(1,end) = 1;
+
+ [X,Y] = meshgrid(x,y);
+
+ figure
+ streamslice(X,Y,u_plot,v_plot,5);
+ title('Velocity Streamlines')
+ ylabel('$\frac{y}{L}$','interpreter','latex')
+ xlabel('$\frac{x}{L}$','interpreter','latex')
+ set(gca,'ylim',[0 1]);
+ axis tight
+
+ u_y_line = u_plot(:,ceil((Npx+2)/2));
+ v_x_line = v_plot(ceil((Npx+2)/2),:);
+
+ format short
+end