帮忙看看问题在哪？谢谢！！！

鸾凤阁

主程序
function mapro

global m
clear,
clf;

%%%%%%%%%%%%%%  参数设置  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k0=3.3*10^7;
g=6.0*10^7;
q=10*10^-3;                                           %tao
p=1.0*10^-14;                                          %ibxt
fs=1.6384*10^6;
ts=0;
tf=0.04;
Ip=0.002;
N=2^16;
h=[k0,g,Ip,p];
%%%%%%%%%%%%%%%% 调用子程序求初态 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

[E0,D0]=initiality(h);
X1=[E0(1);D0(1)];
X2=[E0(2);D0(2)];
X3=[E0(3);D0(3)];
%%%%%%%%%%%%%%%%%%% 调用子程序 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
t=linspace(ts,tf,N-1);
m=input('type into the modulation index m=');
sol=ode45(@vdp2,t,X1);
%%%%%%%%%%%%%%%%%%%%% 求功率谱 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
El=deval(sol,t,1);
D=deval(sol,t,2);
Il=abs(El).^2;
Fl=fftshift(fft(Il));
mFl=abs(Fl(1:N));
f=(0:N-1)*fs/N;
%%%%%%%%%%%%%% 画图 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1)
plot(t,El,'r');
xlabel('time/s');
ylabel('the laser field');
title('时序图')
figure(2)
plot(D,El,'c*');                              
xlabel('population inversion');                                 
ylabel('laser field');
title('相图')
                                                   
figure(3)
plot(t,Il)
xlabel('time/s');
ylabel('the output power/w');
title('激光输出')

figure(4)
plot(f,mFl,'k')
xlabel('frequency in Hz');
ylabel('magnitude of Fl');
title('功率谱')

save data t EL Il D Fl f

子程序1
%calculate initiality
function[x,y]=initiality(h)
x=roots([h(1),-h(4),h(1)+h(1)*h(3)-h(2)*h(3)+h(2),-(1+h(3))*h(4)]);
y=h(1)/h(2)-h(4)/h(2)/x;
子程序2
function subpro=vdp2(t,X)
global  m
subpro=zeros(2,1);
El=X(1);D=X(2);
%%%%%%%%%%%%%%%%%  参数设置  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k0=3.3*10^7;                                %损耗系数
g=6.0*10^7;                                 %增益系数
w=3.5*10^5;                                 %调制角频率
q=10*10^-3;                                 %自发辐射项
p=1.0*10^-14;                               %激光上能级寿命
Ip= 0.002;                                    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
k=k0*(1+m*sin(w.*t));                  %正弦形式的损耗调制

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

subpro(1)=-k*El+g*El.*D+p;
subpro(2)=-((1+Ip+El.^2).^D-Ip+1)/q;
调试时，老出现
???  In an assignment  A(I) = B, the number of elements in B and
I must be the same.
Error in ==> vdp2 at 21
subpro(1)=-k*El+g*El.*D+p;
和
Error in ==> mapro at 33
sol=ode45(@vdp2,t,X1);

ChaChing

Ref: 13F
常见的程序出错问题整理 (eight)
http://forum.vibunion.com/forum/thread-46001-1-1.html