2D convection central difference scheme

here,

Density = 1

Diffusion coefficient=1

When the constant a=10, b=2 in the formula

Steady-state 2D central difference schemeHow should I code?

% Mid_term central difference scheme 
clear all
clc

% global 변수 선언
global AN AS AE AW AP T SC 
NI=4;   % 세로를 4개의 셀로 나눔
NJ=4;   % 가로를 4개의 셀로 나눔
DX=1./(NI-1);
DY=1./(NJ-1);
rho=1.0;	% 밀도를 1.0으로 설정
k=1.0;	   % 확산계수를 1.0으로 설정
u=1.0;   % x방향속도를 1.0으로 설정
v=4.0;   % y방향속도를 4.0으로 설정

% Initialization
X=zeros(NI,NJ);
Y=zeros(NI,NJ);
AP=zeros(NI,NJ);
AN=zeros(NI,NJ);
AS=zeros(NI,NJ);
AE=zeros(NI,NJ);
AW=zeros(NI,NJ);
SC=10-2*(T);  % convection team을 a-b*T 항으로 계산
T=0.*ones(NI,NJ); %모든 내부의 온도 초기값을 0도로 가정

% boundary conditions
%J=1 & NJ
for I=1:NI
T(I,1)=100.;%좌측변 온도를 100도로 설정
T(I,NJ)=0.;%우측변 온도를 0도로 설정
end

%I=1 & NI
for J=1:NJ
T(1,J)=0.;%상변의 온도를 0도로 설정
T(NI,J)=100.;%하변의 온도를 100도로 설정
end

%calculate coefficients
for I=2:NI-1
for J=2:NJ-1
AE(I,J)=-1+0.5*(rho*u);
AW(I,J)=-1-0.5*(rho*u);
AS(I,J)=-1+0.5*(rho*v);
AN(I,J)=-1-0.5*(rho*v);
AP(I,J)=AN(I,J)+AS(I,J)+AE(I,J)+AW(I,J)+(rho*u-rho*u+rho*v-rho*v);
end
end

% call TDMA 
for k=1:1000000 % 100만회 반복계산
  T = TDMA(2,2,NI,NJ,T);   
    if k ~=1
        RESORT=abs(sum(sum(T-T_old)))/(NI*NJ);
        fprintf('ITER= %d RESORT= %e \n', k, RESORT)
    end    
    if k > 1 && RESORT <1.e-6  % 편차가 10-6승 이하면 중지      
        break;
    end
    T_old=T;
end

%-------- Post processing
% grid
x = linspace(0,1,NI);
y = linspace(1,0,NJ); % inherent notation problem 
[X,Y] = meshgrid(x,y);
figure, contourf(X,Y,T,10), colorbar %10+1개의 온도분포선 그림
colormap(gca, jet) % rainbow color map

% 2_d TDMA
function PHI = TDMA(ISTART,JSTART,NI,NJ,PHI)
%	global  AP AN AS AE AW SC ;
     NIM1=NI-1;
     NJM1=NJ-1;
	JSTM1=JSTART-1;
	A(JSTM1)=0.0;
%-----COMMENCE W-E SWEEP
      for I=ISTART:NIM1
          	   C(JSTM1)=PHI(I,JSTM1);
%-----COMMENCE S-N TRAVERCE
	  for J=JSTART:NJM1
%-----ASSEMBLE TDMA COEFFICIENTS
		A(J)=AN(I,J);
		B(J)=AS(I,J);
		C(J)=AE(I,J)*PHI(I+1,J)+AW(I,J)*PHI(I-1,J)+SC(I,J);
		D(J)=AP(I,J);
%%-----CALCULATE COEFFICIENTS OF RECURRENCE FORMULA
	TERM=1./(D(J)-B(J)*A(J-1));
	A(J)=A(J)*TERM;
	C(J)=(C(J)+B(J)*C(J-1))*TERM;
      end
%%-----OBTAIN NEW PHI'S BY BACK SUBSTITUTION
	  for JJ=JSTART:NJM1
	   J=NJ+JSTM1-JJ;
	   PHI(I,J)=A(J)*PHI(I,J+1)+C(J);
      end
      end
end

[YACOUB]

fanction TDMA PHI WERE