format long
%Cavity
%parameters
a=2;
Ra=100;
%discretization and initialization%%%%%%%%%%
N=51; %number of nodes in x- direction
M=(N-1)*a+1;%%number of nodes in y- direction
h=1/(N-1);
k=1/(M-1);
u=zeros(N,M);u_new=u;
T=zeros(N,M);T_new=T;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%coefficients%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

t1=2/(h*h);
t2=2/(k*k)*a*a;
c1=t1/(t1+t2)/2;
c2=t2/(t1+t2)/2;
c4=Ra/2/h/(t1+t2);
c5=k*k/2/(a*a*h*h+k*k);
c6=a*a*h*h/2/(a*a*h*h+k*k);
c7=-a*h*k/8/(a*a*h*h+k*k);
c8=h*h*k*k/2/(a*a*h*h+k*k);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


nr_it=0;    
stop=0;

while (stop~=1)
    nr_it=nr_it+1;
    err_u=0;err_T=0;  
    for i=2:N-1
        for j=2:M-1
             u_new(i,j)=c1*(u(i+1,j)+u(i-1,j))+c2*(u(i,j+1)+u(i,j-1))+c4*(T(i+1,j)-T(i-1,j));
             T_new(i,j)=c5*(T(i+1,j)+T(i-1,j))+c6*(T(i,j+1)+T(i,j-1))+c7*((u(i,j+1)-u(i,j-1))*(T(i+1,j)-T(i-1,j))-(u(i+1,j)-u(i-1,j))*(T(i,j+1)-T(i,j-1)))+c8;
            if abs(u(i,j)-u_new(i,j))>err_u 
                err_u=abs(u(i,j)-u_new(i,j));
            end
            if abs(T(i,j)-T_new(i,j))>err_T 
                err_T=abs(T(i,j)-T_new(i,j));
            end
        
            u(i,j)=u_new(i,j);
            T(i,j)=T_new(i,j);
        end
    end;
   
    T_new(:,1)=T_new(:,2);T_new(:,M)=T_new(:,M-1); % adiabatic condition;
       
    if (err_u<1e-9 & err_T<1e-9)
            stop=1;
        end
        
        
    u=u_new;
    T=T_new;
   if mod(nr_it,500)==0
    fprintf('nr_it=%g err_u=%g err_T=%g\n', nr_it, err_u, err_T);
   end
end
  nr_it
  max(max(T))
  max(max(u))
  x=0:h:(N-1)*h;y=0:k:(M-1)*k;

 figure(1);
  contour(x,y,u',40)
  axis equal
  figure(2);
  contour(x,y,T',40)
  axis equal
% 