Good at Data traning and the various plots
function [psi] = GaussSolver(maxIter,psi,Imax,dE,dN,Jmax,iLe,iTe,jAirfoil,psi_airfoil,err,x_upper,y_upper,x_lower,y_lower,NACA,Alpha)
%Gauss-siedel solver
Omega = 1;global XX YY
cutta = 0;
X=XX;
Y=YY;
psi_old = psi;
for n = 2:maxIter % Iteration Level
psi = psi_old;
for i = 2:Imax-1
for j = 2:Jmax-1
% constants Calculations
%------------------------------------------------------------------------%
r = dE/dN;%abs((X(2,1)-X(1,1))/(Y(1,2)-Y(1,1))); %
if j == jAirfoil(1)
%-------------%
X_i1j2 = (X(i+1,j+2)+X(i+1,j))/2;
X_i_1j2 = (X(i-1,j+2)+X(i-1,j))/2;
Y_i1j2 = (Y(i+1,j+2)+Y(i+1,j))/2;
Y_i_1j2 = (Y(i-1,j+2)+Y(i-1,j))/2;
X_i1j_2 = (X(i+1,j)+X(i+1,j-1))/2;
X_i_1j_2 = (X(i-1,j)+X(i-1,j-1))/2;
Y_i1j_2 = (Y(i+1,j)+Y(i+1,j-1))/2;
Y_i_1j_2 = (Y(i-1,j)+Y(i-1,j-1))/2;
%-------------%
dx_dE_ip = (X(i+1,j)-X(i,j))/dE;
dy_dE_ip = (Y(i+1,j)-Y(i,j))/dE;
dx_dE_im = (X(i,j)-X(i-1,j))/dE;
dy_dE_im = (Y(i,j)-Y(i-1,j))/dE;
dx_dE_jp = (X_i1j2-X_i_1j2)/2/dE;
dy_dE_jp = (Y_i1j2-Y_i_1j2)/2/dE;
dx_dE_jm = (X_i1j_2-X_i_1j_2)/2/dE;
dy_dE_jm = (Y_i1j_2-Y_i_1j_2)/2/dE;