Matrix dimensions must agree,not rendering mesh

Matrix dimensions must agree,not rendering mesh
m=50;
n=50;
t=0.01;
l=0.1;
h=0.1;
r=1000;
w=0.5;
s=0.001;
dumax=0.001;
dvmax=0.001;
dpmax=0.001;
dx=l/(m+1);
dy=l/(n+1);
ou=zeros(m);
ou(1,:)=1;
ov=zeros(m);
op=zeros(m);
op(:,:)=10315;
x=zeros(m,1);
y=zeros(1,n);
for i=1:m-2
x(i)=(i-1)*dx;
end
for j=1:n-2
y(j)=(j-1)*dy;
end
for i=2:m-1
for j=2:n-1
oA1(i,j)=(ou(i+1,j)-2*ou(i,j)+ou(i-1,j))/dx^2;
oA2(i,j)=(ou(i,j+1)-2*ou(i,j)+ou(i,j-1))*dx^2;
oA3(i,j)=(ou(i+1,j)^2-ou(i,j)^2)/dx;
oA4(i,j)=(ou(i,j+1)*ov(i,j+1)-ou(i,j)*ov(i,j))/dy;
oA(i,j)=1/r*(oA1(i,j)+oA2(i,j))-oA3(i,j)-oA4(i,j);
su(i,j)=t*(oA(i,j)-(t/dx)*(op(i+1,j)-op(i,j))+ou(i,j));
end
end
for i=2:m-2
for j=2:n-2
oB1(i,j)=(ov(i+1,j)-2*ov(i,j)+ou(i-1,j))/dx^2;
oB2(i,j)=(ov(i,j+1)-2*ov(i,j)+ou(i,j-1))*dx^2;
oB3(i,j)=(ov(i+1,j)^2-ov(i,j)^2)/dx;
oB4(i,j)=(ov(i,j+1)*ou(i,j+1)-ov(i,j)*ou(i,j))/dy;
oB(i,j)=1/r*(oB1(i,j)+oB2(i,j))-oB3(i,j)-oB4(i,j);
sv(i,j)=t*(oB(i,j)-(t/dx)*(op(i,j+1)-op(i,j))+ov(i,j));
end
end
for i=2:m-1
for j=2:n-1
a=2*t*(1/dx^2+1/dy^2);
b=t/dx^2;
c=t/dy^2;
d(i,j)=1/dx*(ou(i,j)-ou(i-1,j))+1/dy*(ov(i,j)-ov(i,j-1));
p1(i,j)=(b*op(i+1,j)+op(i-1,j)+c*op(i,j+1)+c*op(i,j-1)-d(i,j))/a;
p1(i,j)=w*p1(i,j)+(1-w)*p1(i,j);
end
end
for i=1:m-2
for j=1:n-2
nu(i,j)=ou(i,j)-t*(p1(i+1,j)-p1(i,j))/dx;
nv(i,j)=ov(i,j)-t*(p1(i,j+1)-p1(i,j))/dy;
np(i,j)=op(i,j)+p1(i,j);
end
end
for i=1:m-2
for j=1:n-2
du=abs(nu(i,j)-ou(i,j));
dv=abs(nv(i,j)-ov(i,j));
dp=abs(np(i,j)-op(i,j));
if du>dumax&dv>dvmax&dp>dpmax;
ou=nu;
ov=nv;
op=np;
return
else
break
end
end
end
高手可不可以帮我把这个程序实现循环啊,我弄了好长时间了,就是不循环!

load DATA.DAT; N=length(DATA) x=DATA(1:N,1); y=DATA(1:N,2); z=DATA(1:N,3); minx = min(x); maxx = max(x); miny = min(y)