Matlab Code

clc clear all %% M=1000; mMin=1; mMax=M; N=100; nMin=1; nMax=N; xMin=0; xMax=0.3; deltaX=(xMax-xMin)/(M-1); %X-Parameter yMin=0; yMax=0.0047; deltaY=(yMax-yMin)/(N-1); %Y-Parameter Ux=1; %FreeStream Velocity (m/s) mu=10^-6; %Kinematic Viscosity (m^2/s) rho=1000; %Density (kg/m^3) x=0.01; %% for i= mMin:mMax %Define Domain u and v from Boundary Conditions and Initial Conditions for j= nMin:nMax u(nMin,i)=0; v(nMin,i)=0; u(nMax,i)=Ux; v(nMax,i)=0; u(j,mMin)=Ux; v(j,mMin)=0; alpha(j,i)= (mu*deltaX)/((u(j,i))*deltaY^2); alpha2(j,i)=1+2*alpha(j,i); beta(j,i)=(v(j,1)*deltaX)/(u(j,i)*(2*deltaY)); end end %% for j=nMin+1:nMax-2 A=full(gallery('tridiag',N-2,-alpha(j,1),alpha2(j,1),-alpha(j,1))); %Tridiagonal Matrix A*B=X B(j,1)=u(j,1)-beta(j,1)*(u(j+1,1)-u(j-1,1)); end %% % Tridiagonal Matrix Solver X2 = zeros(length(B),1); c = B; d = diag(A); b = zeros(nMax-2,1); b(2:end) = diag(A,-1); a = zeros(nMax-2,1); a(1:end-1) = diag(A,+1); test = zeros(nMax-2,1); test(2:end) = b(2:end)./d(1:end-1); d(2:end) = d(2:end) - test(2:end).*a(1:end-1); c(2:end) = c(2:end) - test(2:end).*c(1:end-1); X2(end) = c(end)/d(end); for z = N-3:-1:1 X2(z) = (c(z) - a(z)*X2(z+1))/d(z); end %% X1=0; %n=1 values X3=Ux; %n=N values u(:,mMin+1)=[X1;X2;X3]; %Update Domain u %% %Update Domain v for j=nMin+1:nMax-1 v(j,2)=v(j-1,2)-(deltaY/(2*deltaX))*(u(j,2)-u(j,1)+u(j-1,2)-u(j-1,1)); end v(nMin,2)=0; v(nMax,2)=0;