clear?
clc
RE=5000;
niu=1.139*10^（-6）;
TN=500000;
err=0.0001;
rho=0.05;
L=1;
h=0.01;
N=L/h+1;
u_0=niu*RE/L;
Tw=5000;
E1=0;
E2=1;
counter=0;
U=zeros（NN）;%流量
U0=zeros（NN）;
Z=zeros（NN）;%渦量
Z0=zeros（NN）;
Vx=zeros（NN）;
Vy=zeros（NN）;

%t=0時刻初始值即U?Z?VX?VY均為0;
for?t=1:TN
????%求邊界上的渦量值
????%DC
????for?i=2:N-1
????????Z（iN）=2*（（U（iN）-U（iN-1）-u_0*h）/h^2）;
????end
????%AB
????for?i=2:N-1
????????Z（i1）=2*（（U（i1）-U（i2））/h^2）;
????end
????%AD
????for?j=2:N-1
????????Z（1j）=2*（（U（1j）-U（2j））/h^2）;
????end
????%BC
????for?j=2:N-1
????????Z（Nj）=2*（（U（Nj）-U（N-1j））/h^2）;
????end
????%四個頂點的渦量值
????Z（11）=0.5*（Z（12）+Z（21））;
????Z（1N）=0.5*（Z（1N-1）+Z（2N））;
????Z（N1）=0.5*（Z（N-11）+Z（N2））;
????Z（NN）=0.5*（Z（NN-1）+Z（N-1N））;
????
????%由2.7.24求內點的新渦量值
????for?i=2:N-1
????????for?j=2:N-1
????????????u=max（max（abs（Vx）））;
????????????v=max（max（abs（Vy）））;
????????????dt=min（h^2/4/niuh/（u+v））;
????????????Z（ij）=Z（ij）+（Vx（ij）*（Z（i-1j）-Z（i+1j））+Vy（ij）*（Z（ij-1）-Z（ij+1）））...
????????????????*（dt/2/h）+（niu*dt/h^2）*（（Z（i+1j）-2*Z（ij）+Z（i-1j））+（Z（ij+1）-2*Z（ij）