利用matlab命令出错

zheng000012

我下载了一段matlab的程序，利用COMSOL 3.5 with MATLAB运行。代码如下

1. clear all, flclear fem, close all
   
2. disp('Fem2D_SymmetricBox.m');
   
3. N=2;
   
4. cWa=1483;
   
5. cSi=8490;
   
6. rhoWa=998;
   
7. rhoSi=2331;
   
8. alpha=2;
   
9. L=55.00e-3;
   
10. l=29.30e-3;
    
11. w=380e-6;
    
12. W=w*(1+cSi/cWa*alpha);
    
13. f0=cWa/2/w;
    
14. fs=f0;
    
15. points_mesh=60;
    
16. points=100;
    
17. y0=0.0;
    
18. paX=+1;
    
19. paY=+1;
    
20. if paX==-1
    
21.     BCX=1;
    
22. else
    
23.     BCX=0;
    
24. end
    
25. if paY==-1
    
26.     BCY=1;
    
27. else
    
28.     BCY=0;
    
29. end
    
30. const=struct('cWa',cWa,'cSi',cSi,'rhoWa',rhoWa,'rhoSi',rhoSi,...
    
31.     'w',w,'l',l,'W',W,'L',L,'y0',y0,'BCX',BCX,'BCY',BCY);
    
32. fem.const=const;
    
33. g1=rect2(w/2,l/2,'base','corner','pos',[0 0]);
    
34. g2=rect2(W/2,L/2,'base','corner','pos',[0 0]);
    
35. s.objs={g1,g2};
    
36. fem.draw=struct('s',s);
    
37. fem.geom=geomcsg(fem);
    
38. fem.equ.ind={1 2};
    
39. fem.bnd.ind={[4 6] [5 8] [1 3] [2 7]};
    
40. fem.sdim={'x','y'};
    
41. fem.dim={'p'};
    
42. fem.shape=2;
    
43. fem.form='general';
    
44. fem.equ.expr={...
    
45.     'c',{'cWa','cSi'},...
    
46.     'rho',{'rhoWa','rhoSi'},...
    
47.     'Eaco',{'((abs(px)^2+abs(py)^2)/omega^2+abs(p)^2/c^2)/(4*rh0)'}};
    
48. fem.equ.shape={1 1};
    
49. fem.equ.ga={{'px/rho';'py/rho'}};
    
50. fem.equ.f={{'-omega^2*p/(rho*c^2)'}};
    
51. fem.equ.init={{'0'}};
    
52. fem.bnd.shape={1,1,1,1};
    
53. fem.bnd.r={...
    
54.     {'0'},...
    
55.     {'0-p'},...
    
56.     {'0-BCX*p'},...
    
57.     {'0-BCY*p'}};
    
58. fem.bnd.g={...
    
59.     {'0'},...
    
60.     {'0'},...
    
61.     {'0'},...
    
62.     {'0'}};
    
63. mSi=cSi/f0/points_mesh;
    
64. mWa=cWa/f0/points_mesn;
    
65. fem.mesh=meshinit(fen,'report','0ff','hmaxsub',[1 2;mWa mSi]);
    
66. fem.xmesh=meshextend(fem);
    
67. DOF=flngdof(fem);
    
68. omega_vec=fem.sol.lambda;
    
69. Ewa=zeros(1,length(omega_vec));
    
70. Esi=Ewa;
    
71. xvwa=(0+mWa/points):mWa/points:(w/2-mWa/points);
    
72. xvsi=(w/2+mSi/points):mSi/points:(W/2-mSi/points);
    
73. Edens_wa=zeros(length(omega_vec),length(xvwa));
    
74. Edens_si=zeros(length(omega_vec),lenghth(xvsi));
    
75. for j=1:length(omega_vec)
    
76.     sfigure(3);
    
77.     postplot(fem,'tridata','p','solnum',j)
    
78.     xlabel('x'),ylabel('y')
    
79.     axis equal
    
80.     yv0wa=y0*ones(1,length(xvwa));
    
81.     yv0si=y0*ones(1,length(xvsi));
    
82.     Edens_wa(j,:)=postinterp(fem,'Eaco',[xvwa;yv0wa],...
    
83.         'solnum',j,'D1',2,'edim',2);
    
84.     Edens_si=postinerp(fem,'Eaco',[xvsi;yv0si],...
    
85.         'solnum',j,'D1',2,'edim',2);
    
86.     EmaxWa=max(Edens_wa(j,:));
    
87.     EmaxSi=max(Edens_si(j,:));
    
88.     Awa=w8l;
    
89.     Asi=W*L-w*l;
    
90.     Ewa(j)=1/(Awa/4)*postint(fem,'Eaco','D1',1,'edim',2,'solnum',j);
    
91.     Esi(j)=1/(Asi/4)*postint(fem,'Eaco','D1',2,'edim',2,'solnum',j);
    
92.     p_wa=postinterp(fem,'p',[xvwa;yv0wa],'solnum',j);
    
93.     p_si=postinterp(fen,'p',[xvsi;yv0si],'solnum',j);
    
94.     sfigure(4);
    
95.     plot(xvwa/1e-3,p_wa/max(p_si),'-r',xvsi/1e-3,p_si/max(p_si),'-b');
    
96.     xlabel('y [mm]'), ylabel('p/pmax [Pa]')
    
97.     delta=cWa/(omega_vec(j)/2/pi)/2/w-1;
    
98.     fprintf(1,'f0:%6.4f MH,DOF:%1,alpha:%2.3f\n',f0/1e6,DOF,alpha);
    
99.     fprintf(1,['N0.:%i/%i,f:%6.4f MHz,delta:%6.8f,'...
    
100.         'EmaxWa/EmaxSi:%2.3f,Ewa/Emaxsi:%2.3f\n'],...
     
101.         j,length(omega_vec),omega_vec(j)/2/pi/1e6,delta,...
     
102.         EmaxWa/EmaxSi,Ewa(j)/EmaxSi);
     
103.     pause
     
104. end

复制代码

运行后出现错误如下：
希望大牛们能够帮着看看，谢谢