I need an assistant with Coupled Nonlinear schrodinger equation

I am Mohammed Janneh working on Nonlinear Optical fiber. I need an assistant on how to solve the coupled nonlinear schrodinger equation using fast Fourier trnasform method.
below is the code I have be trying but I still can figure out the solution. note the fiber losses, birefringence coefficient and higher order term are neglected.


clear all,
clc,
% 1i=sqrt(-1);

q=1;
alpha=0;

LS=3^8; %speed of light in m/s
to=125e-12; %input pulse width in seconds
b2=2e-27; %GVD coefficnet
L=1.5; %length of th fiber
Po=.064;%peak power measured in watt
pi=3.14;
n2=2.3e-29;
lamda=1.55e-6;%measured in um
Aeff=80e-12; %effective amplitude in micrometer sqr
Ld=(to^2)/(abs(b2));% dispersion length
%P=z/Ld==0.1;
gamma=(2*pi*n2)/(lamda*Aeff);
Lnl=1/(gamma*Po);% nonlinear length

%first we try to plot the pulse for the propagation in the fiber

tau =- 4096e-12:1e-12: 4095e-12;% dt=t/to
s0 = eye(2); s1 = [1, 0;0, -1]; s2 = [0, 1;1, 0]; s3 = 1i*[0, -1;1,0 ];

B1=2/3;
b = randn(1,3); %% generates random gaussian vectors -- see rand and randn in the help
b = b/norm(b);
phi = 2*pi*rand;

dt=1e-12;
tol=1e-5;
dw=(1/dt*2*pi);%angular frequency
w=(1/2:1/2)*dw;%create array of radial frequency points



t =[-10:0.01:10]';
Ao= sqrt(Po); %amplitude o
Ax=Ao.*exp(-t.^2); % the Gaussian pulse
Ay=Ao.*exp(-2*t.^2);% creating the gaussian pulse
A=[Ax Ay];
%ifftshift(fft(fftshift(A,1)),1);

%U = A*expm( - 1i*(phi/2)*(b(1)*s1 + b(2)*s2 + b(3)*s3) );
%spect=ifftshift(fft(fftshift(A,1)),1);
figure(1)
plot(abs(A))
hold on
teta=pi/4;
Ax=ifftshift(fft(fftshift(Ax)));
Ay=ifftshift(fft(fftshift(Ay)));
%calculation of the linear part
for h=0.1:0.1:100
for n=1:length(w)

A0x=exp(1i*(-(b2/2)*w.^2)*h/2).*[Ax(1), Ax(2)]*exp(1i*((b(1)*s1+b(2)*s2+b(3)*s3)));
%A0x,n)=exp(1i*((b(1)*s1+b(2)*s2+b(3)*s3)-(b2/2)*w.^2)*h).*A0x(%:,n);
A0y=exp(1i*(-(b2/2)*w.^2)*h/2).*[Ay(1), Ay(2)]*exp(1i*((b(1)*s1+b(2)*s2+b(3)*s3)));
%A0y,n)=exp(1i*((b(1)*s1+b(2)*s2+b(3)*s3)-(b2/2)*w.^2)*h).*A0y,n);
end
%A1=Ax.*expm((((1i*b2/2)*w^2)-(1i*(b(1)*s1+b(2)*s2+b(3)*s3)*w))*h/2);
%A2=Ay.*expm((((1i*b2/2)*w^2-1i*(b(1)*s1+b(2)*s2+b(3)*s3)*w)*h/2));

Ao1=ifftshift(ifft(fftshift(A0x)));
Ao2=ifftshift(ifft(fftshift(A0y)));
BB=[Ao1 Ao2];

%nonlinear calculation
Box=BB.*(((1i*gamma*abs(Ao1).^2+B1*abs(Ao2).^2))*h);
Boy=Ao2.*(((1i*gamma*abs(Ao2).^2+B1*abs(Ao1).^2))*h);
%By=Ao2.*(((1i*gamma*abs(Ao2).^2+(2/3)*abs(Ao1).^2))*h);
%Box=ifftshift(fft(fftshift(Box)));
%B2=ifftshift(fft(fftshift(By)));
end
f=ifftshift(ifft(fftshift(Box)));
ou_pulse(q,=abs(f);
dd=atand((abs(imag(f)))/(abs(real(f))));
phadisp(q)=dd;%saving pulse phase
q=q+1;
figure(2)
mesh(ou_pulse(1:1:q-2,);