clear

     NN =input(  'NN=  ');
     beta=6; p=20; N=2^NN; lamda=0.065; 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;

     u=[2.*sech(xx).^2];

     time=dt;

     while time < T 

            K1=codeKdV_RHS(u,N,beta,mu,p);
            K2=codeKdV_RHS(u+0.5*dt*K1,N,beta,mu,p);
            K3=codeKdV_RHS(u+0.5*dt*K2,N,beta,mu,p);
            K4=codeKdV_RHS(u+dt*K3,N,beta,mu,p); 

          u=u+(dt/6)*(K1+2*K2+2*K3+K4); 

          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, '-')