How to code split-step Fourier method in Matlab?

Question

I need to solve the following split-step equation in Matlab:

$$u(x+\Delta x, z)=\text{exp}\bigg(ik(n^2-1)\frac{\Delta x}{2}\bigg)F^{-1}\bigg(\text{exp}\bigg[-ip^2\frac{\Delta x}{2k}\bigg]F\{u(x,z\}\bigg)$$

I have an initial field $u(0,z)$, which is a vector, and I want the Matlab to find $u(x+\Delta x, z)$ for some range $x$. Effectively i want something like this:

In my head when I code something like this:

xmax=2600; 
delx=50;  % step size in x
x = delx:delx:xmax;

 for i = delx:delx:xmax;   
   c1=fftshift(fft(u0z));
   c2=exp(-1j*(p^2.*x/(2*k))).*c1;
   c3=ifft(ifftshift(c2));
   u0z=exp(1j*k.*(n.^2-1).*x./2).*c3;
 end

The program should use my initial vector $u0z$, feed it and find the value of $u(x+\Delta x, z)$. Then it should feed the new value and use it to find $u(x+2\Delta x, z)$ and so on, but it doesn't. It just gives me a $1 x n$ row vector at the end.

I am sure there is something wrong with how I am implementing code, and perhaps even how I use the fft commands, but I don't really know Matlab. Would anyone be willing to help? Thanks.

P.S. I could paste my whole code if required, but it is long and messy.

EDIT:

for i=1:1:length(x)
   c1=fftshift(fft(u0z));
   c2=exp(-1j*(p^2*x(i)/(2*lamda))).*c1;
   c3=ifft(ifftshift(c2));
   u0z=exp(1j*lamda*(n(i)^2-1)*x(i)/2).*c3;
end