format long g;
tic
N=51;
h=1/(N-1);
%tol=1e-8;
xplot=0:h:1;
yplot=0:h:1;

Re=10; 

u=zeros(N,N);O=zeros(N,N);
uo=zeros(N,N);Oo=zeros(N,N);
Utop=-1;

Oo(:,N)=2*(u(:,N)-u(:,N-1))/h/h-2*Utop/h;

stop=1;
nr_it=0;
rO=0.9;
ru=1;
while stop==1
      nr_it=nr_it+1;
   
      for i=2:N-1
       for j=2:N-1
           
        O(i,j)=0.0625*Re*((Oo(i,j+1)-O(i,j-1))*(uo(i+1,j)-u(i-1,j))-(Oo(i+1,j)-O(i-1,j))*(uo(i,j+1)-u(i,j-1)))+...
            0.25*(Oo(i+1,j)+O(i-1,j)+Oo(i,j+1)+O(i,j-1));
        O(i,j)=rO*O(i,j)+(1-rO)*Oo(i,j);
     
        u(i,j)=0.25*(uo(i+1,j)+u(i-1,j)+uo(i,j+1)+u(i,j-1))+0.25*h*h*O(i,j);
       % u(i,j)=ru*u(i,j)+(1-ru)*uo(i,j);
       end;
      end;

     O(1,:)=-2*(u(2,:)-u(1,:))/h/h;
     O(N,:)=2*(u(N,:)-u(N-1,:))/h/h;
     O(:,1)=-2*(u(:,2)-u(:,1))/h/h;
     O(:,N)=2*(u(:,N)-u(:,N-1))/h/h-2*Utop/h;
     
    errU=norm(u-uo);errO=norm(O-Oo);
      
      if mod(nr_it,1000)==0
         fprintf('nr_it=%d  errU=%g  errO=%g\n', nr_it, errU, errO);
         if (errU<1e-7)&(errO<1e-7)
             stop=0;
         end
        
      end
      uo=u;Oo=O;
        
end;

 contour(xplot,yplot,u',100);
axis equal
max(max(abs(u)))
toc