clear; close 

     NN =input(  'NN=  ');
     beta=6; p=20; N=2^NN; lamda=0.0323; mu=1; T=1.0;

     dx=2*p/N; dt=lamda*(dx)^3;

     j=[0:1:N-1];  x=[2*pi*j/N];  xx=(x/pi-1)*p;

     u0=[2.*sech(xx).^2];
     u1=[2.*sech(xx-4*dt).^2];

     time=dt;

     while time < T 
        rhs=codeKdV_RHS(u1,N,beta,mu,p);
        u=u0+2*dt*rhs;
        u0=u1;  u1=u;
        time=time+dt;
     end

     figure

     v=u'

     uexact=[2.*sech(xx-4*(time-dt)).^2];

     err=max(abs(uexact-u))

     fprintf(1,      '%9.0f          %10.4e  \n', [N;err]);

     axis([-5, 10, 0, 2])
     plot (xx, real(u), 'x', xx, uexact, '-')