L = (pi/2)*t0^2/abs(beta);
h = L/1000;
gama= 3;
E0 = (1/t0)*sqrt(abs(beta)/gama);
tau=(tf-ti)/(K-1);
alpha=1i*beta*h/tau^2;
A=zeros(K-2,K-2);
% Compute A
for i=1:K-2
    for j=1:K-2
        if i==j-1   ,A(i,j)= -1*alpha;  end
        if i==j     ,A(i,j)= 2*alpha;   end
        if i==j+1   ,A(i,j)= -1*alpha;  end
    end
end
t=(ti+tau:tau:tf-tau)';
u0=N*E0*sech(t/t0);
v0=1i*2*gama*h*((conj(u0)).*u0).*u0;
u1= 0.5*A*u0 + 0.5*v0 + u0;
u=zeros(K,(L/h));
u(2:K-1,1)=u0;
u(2:K-1,2)=u1;
j=1;
period=1;
while j <= period*(L/h)-2
    j=j+1;
    v=1i*2*gama*h*(conj(u(2:K-1,j)).*u(2:K-1,j)).*u(2:K-1,j);
    u(2:K-1,j+1)=A*u(2:K-1,j) + v + u(2:K-1,j-1);
end
U2=((conj(u)).*u)';

clf
[r c]=size(U2);
x=linspace(ti/t0,tf/t0,c);
y=linspace(0,period,r);
mesh(x,y,U2);
view(30,30)
axis tight
shading interp
xlabel('Time ( t / t_0 )')
ylabel('Distance ( L / L_0 )')
title('Intensity')

--
Reference

S. Lopez-Aguayo, J.P. Ochoa-Ricoux, and J.C. Gutierrez-Vega "Exploring the behavior of solitons on a desktop personal computer" REVISTA MEXICANA DE FISICA E 52 (1) 28–36.