3rd order nonlinear PDE

if true
  % clear all; close all; m1=1;m2=0.1;m3=0.1;xL=-10; xR=10; Nx=50; t0=0; tf=0.4; Nt=25;
h=(xR-xL)/Nx; 
k=(tf-t0)/Nt; 
r=m2*k/h^2;
s=m3*k/2*h^3;
x=xL:h:xR; 
u=zeros(Nx+1,Nt); 
%if exact solution is 
f=(-6*m3^2/25*m2)*(1 + tanh(m3*x/10*m2) + 0.5*sech(m3*x/10*m2)^2);
for j=1:Nt t=j*k;
if j==1
for i=2:Nx u(i,j)=f(i)+r*(f(i+1)-2*f(i)+f(i-1))-s*(f(i+2)-2*f(i+1)+2*f(i-1)-f(i-2))-k/h*f(i)*(f(i+1)-f(i)); end;
u(1,j)=0; u(Nx+1,j)=0;
else
for i=2:Nx
u(i,j)=u(i,j-1)+r*(u(i+1,j-1)-2*u(i,j-1)+u(i-1,j-1))-s*(u(i+2,j-1)-2*u(i+1,j-1)+2*u(i-1,j-1)-u(i-2,j-1))-k/h*u(i,j-1)*(u(i+1,j-1)-u(i,j-1));
end;
u(1,j)=0; u(Nx+1,j)=0;
end
end; T=1*k:k:Nt*k;
figure(1);
surf(x,T,u'); xlabel('x'); ylabel('y'); zlabel('u');
title('Numerical solution');
figure(2);
hold on; plot(x,u(:,1),'k-',x,u(:,10),'k--',x,u(:,Nt),'k:','LineWidth',2);
title('Solution profile at t_k= 0.016, 0.16, 0.4');
legend('t=0.016','t=0.16','t=0.4');
figure(3);
G=plot(x,u(:,1),'LineWidth',3,'erase','xor');
for i=2:Nt set(G,'xdata',x,'ydata',u(:,i)); pause(0.1); end
title('Animation of solution'); xlabel('x'); ylabel('u(x,t_k)');
end

[Bjorn Gustavsson]

Please edit your post so that the code section is done with the CODE-markup.