Double Slit Diffraction Pattern

Suppose that we have two slits, the distance between them is d=10 and their width a=0.1. My initial condition is a Gaussian. In order to compute the problem, our teacher told us to use the split-step Fourier method. My try at the code is shown below. For some reason, I don't get any diffraction patterns and instead only two straight lines at the points of my po function. What can I do?

clc; clear all; format long;
N=512; L=100; x=(L/N)*(-N/2:N/2-1); dx=x(2)-x(1);
k=(2*pi/L)*[0:N/2-1 -N/2:-1]; 
w=10; uin=exp(-(x/w).^2); % initial condition
%figure(1); plot(x,abs(uin),'r');

S=1000; zmax=20; h=zmax/S; z=0; count=0; dim=0; savestep=10; u0=uin; 
for m=1:S
po=rectangularPulse(0,1/10,x) + rectangularPulse(101/10,102/10,x);
u1=ifft(fft(u0).*exp(-1i*0.5*h*k.^2));   % First step of SSF
u2=exp(1i*h*po).*u1;                     % Second step of SSF
u3=ifft(fft(u2).*exp(-1i*0.5*h*k.^2));   % Second step of SSF
u0=u3; z=z+h; count=count+1; 
if (count==savestep)
    dim=dim+1; field(dim,:)=u0; pot(dim,:)=po; z1(dim)=z; count=0;
end;
end; 
[X,Z]=meshgrid(x,z1);
figure(2); surf(X,Z,abs(field)); shading interp; view([0 90]);  axis square; axis tight;
figure(3); surf(X,Z,pot); shading interp; view([0 90]);  axis square; axis tight;

[Sam Chak]

@Apostolis

Where is your < m > in the loop? I guess, once you fix it, you can find a way to generate the surface plot.