色散

function dP=fwm(z,P)
dP=zeros(8,1);
temp1=-6.77e4*z+P(7)+P(8)-P(5)-P(6);
temp2=(P(1)*P(2)*P(3)*P(4))^0.5;
dP(1)=-500*P(1)-4/3*1522000*temp2*sin(temp1);
dP(2)=-500*P(2)-4/3*1522000*temp2*sin(temp1);
dP(3)=-500*P(3)+4/3*1556000*temp2*sin(temp1);
dP(4)=-500*P(4)+4/3*1487000*temp2*sin(temp1); 
dP(5)=2/3*1522000*temp2/P(1)*cos(temp1);
dP(6)=2/3*1522000*temp2/P(2)*cos(temp1);
dP(7)=2/3*1556000*temp2/P(3)*cos(temp1);
dP(8)=2/3*1487000*temp2/P(4)*cos(temp1);
end 

m=1;
y1=zeros(2001,1).';
y2=zeros(2001,1).';
for a=0:0.001:2
P0=7.05;
[Z,P]=ode45('fwm',[0,0.00125],[0.001;0.001;P0;P0;-a*pi;-a*pi;-0.25*pi;-0.25*pi]);
A1=P(:,1).';
A2=P(:,2).';
A3=P(:,3).';
A4=P(:,4).';
A5=P(:,5).';
A6=P(:,6).';
A7=P(:,7).';
A8=P(:,8).';
y1(m)=10*log(A1(1,length(A1))/0.001);
y2(m)=(2.51e3*0.00125+A7(1,l