N=input（‘輸入N：‘）;?
M=input（‘輸入M：‘）;??
K=input（‘輸入K：‘）;??
R=1;????%單元劃分數（2N+M）*K圓柱底面半徑為1；
%劃分結構化網格
L=10*R;
rr=linspace（0RN+1）;
ll=linspace（0LM+1）;
thet=linspace（02*piK+1）;
x=zeros（2*N+M+1K+1）;
y=zeros（2*N+M+1K+1）;
z=zeros（2*N+M+1K+1）;
P=cell（2*N+MK）;
A=cell（2*N+M+1K+1）;
delta_S=zeros（2*N+MK）;
delta_r=R/N;
delta_l=L/M;
delta_thet=2*pi/K;

for?i=1:N;
????x（i:）=0;
????for?j=1:K+1;
????????z（ij）=rr（i）*sin（thet（j））;
????????y（ij）=rr（i）*cos（thet（j））;
????????A{ij}=[x（ij）y（ij）z（ij）];
????end
end
for?i=N+M+1:2*N+M+1;
????x（i:）=L;
????for?j=1:K+1;
????????z（ij）=rr（2*N+M+2-i）*sin（thet（j））;
????????y（ij）=rr（2*N+M+2-i）*cos（thet（j））;
????????A{ij}=[x（ij）y（ij）z（ij）];
????end
end
for?i=N+1:N+M;
????x（i:）=ll（i-N）;
????for?j=1:K+1;
????????z（ij）=R*sin（thet（j））;
????????y（ij）=R*cos（thet（j））;
????????A{ij}=[x（ij）y（ij）z（ij）];
????end
end

for?j=1:K;
?????P{1j}=1/3*（A{1j}+A{2j}+A{2j+1}）;
?????P{2*N+Mj}=1/3*（A{2*N+Mj}+A{2*N+Mj+1}+A{2*N+M+1j}）;
end

for?i=1:N;
????delta_S（i:）=sin（delta_thet）*（rr（i+1）^2-rr（i）^2）/2;
????delta_S（2*N+M-i+1:）=delta_S（i:）;
end
sw=sin（delta_thet/2）*R*2*delta_l;
for?i=N+1:N+M;
????delta_S（i:）=sw;
end
for?i=2:2*N+M-1;
????for?j=1:K;
????????P{ij}=0.25*（A{ij}+A{ij+1}+A{i+1j}+A{i+1j+1}）;
????end
end?
surf（xyzz）
xlabel（‘x‘）
ylabel（‘y‘）
zlabel（‘z‘）
axis?image

PP=2*N+M;
%開始進行單元計算，解方程組
%設遠方速度為1，x正方向
a=zeros（PP*KPP*K）;
b=zeros（PP*K1）;
c=zeros（PP*KPP*K）;
sigma=zeros（PP*K1）;
PHI=zeros（PP*K1）;
r=zeros（PP*KPP*K）;

%求圓柱m_11
n_p=cell（PP*K1）;
PHI_1=zeros（PP*K1）;

for?k=