clf;
n=50;

a=1;
c=2.1; % try changing it to 1.9

Tfinal=13000;
Tstep=10;

x=linspace(0,n-1,n);
LEN=n;

% Initial conditions
u=sin(2*pi/LEN*x')*0.01-sin(6*pi/LEN*x')*0.2+1;
U=u';

tnext=0;
thist=[]; uhist=[];
t=0;
dt=0.1;
while t<Tfinal

   t=t+dt;

   idxp=[2:n,1];
   idxm=[n,1:n-1];
   Tot=U(idxp)+U(idxm);
   etap=U(idxp)./Tot;
   etam=U(idxm)./Tot;

   rhsU=a*(U(idxp)+U(idxm))-(2*a+c)*U+c*U(idxp).*etam(idxp)+c*U(idxm).*etap(idxm);
   U=U+dt*rhsU;

   if t>tnext,
      tnext=tnext+Tstep;
      
      plot(x,U,'ko-'); 
      title(sprintf('t=%g',t));
      thist(end+1)=t;
      uhist(end+1,:)=u;
      drawnow;
   end;
end;

%%
%save('full_discr.mat', 'thist','uhist','x');
