仅供参考
clc;
clear all;
a=1;
b=2;
q0=1/32;
D=1/pi^6;
z=0;
Nx=10;
Ny=10;
xs=linspace(0,1,Nx);
ys=linspace(0,2,Ny);
for i=1:Nx
for j=1:Ny
x=xs(i);
y=ys(j);
for m=1:2:100
for n=1:2:100
fun =16*q0*sin(m*pi*x/a)*sin(n*pi*y/b)./(pi^6*D*m*n*(m^2/a^2+n^2/b^2)^2);
z=z+fun;
end
end
zs(i,j)=z;
end
end
[xs1,ys1]=meshgrid(xs,ys);
surf(xs,ys,zs)
