myalpha=0.25;   %Nonlinear parameter
recurrence=1; % show recurrence or kinks?

if(recurrence)
 N=32*1;         %Number of particles must be a power of 2
 TMAX=10000; DT=10;
else
 N=32*2;         %Number of particles must be a power of 2
 TMAX=100; DT=1;
end
%TMAX=10000; DT=20; % to show recurrence
tspan=[0:DT:TMAX];

options=odeset('Reltol',1e-4,'OutputFcn','odeplot','OutputSel',[1,2,N]); % Test different tolerances, changing Reltol
options=odeset('Reltol',1e-4);

for ii=1:N,
 if(recurrence)
  a=1; b(ii)=a*sin(pi*ii/(N+1)); b(ii+N)=0;  % FPU initial condition
  %a=1; b(ii)=a*sin(pi*N*ii/(N+1)); b(ii+N)=0; % Zabusky-Deem init. cond.
  omegak2(ii)=4*(sin(pi*ii/2/N))^2;   % Mode Frequencies
 else % kinks
  k=0.2; sk=(sinh(k))^2; ek=exp(-k); i1=ii-N/4; i2=i1-N/2; %Solitons
  %b(ii)=0; b(ii+N)=0;
  b(ii)=     -0.5/myalpha*log((1+exp(2*k*(i1-1)))/(1+exp(2*k*i1)));  % Kink
  b(ii)=b(ii)+0.5/myalpha*log((1+exp(2*k*(i2-1)))/(1+exp(2*k*i2)));  % Anti-kink
  b(ii+N)=       -sk*ek/myalpha/cosh(k*i1)/(exp(-k*i1)+exp(k*i1)*exp(-2*k)); % Kink vel
  b(ii+N)=b(ii+N)+sk*ek/myalpha/cosh(k*i2)/(exp(-k*i2)+exp(k*i2)*exp(-2*k)); % Anti-Kink vel
 end
end

IC=b';

[T,Y]=ode45(@(t,y) fpu1(t,y,N,myalpha),tspan,IC,options);  % Time integration

mymovie=1;
if(mymovie)
 figure(1);clf
 mm=min(min(Y(:,1:N)));
 MM=max(max(Y(:,1:N)));
 for ii=1:length(T)
 plot(1:N,Y(ii,1:N),'o-')
 axis([1 N mm MM])
 title(['time=',num2str(T(ii))])
 drawnow
 end
end %if(mymovie)


 figure(2);clf
 surf(1:N,T,Y(:,1:N))
 axis tight
 xlabel('n')
 ylabel('t')
 zlabel('u_n')

myenergy=0;
if(myenergy)
 for iiT=1:(TMAX/DT),
   TiiME(iiT)=iiT*DT*sqrt(omegak2(1))/2/pi;  % Time iteration loop
   YX(iiT,1:N+1)=[0 Y(iiT,1:N)];  
   YV(iiT,1:N+1)=[0 Y(iiT,N+1:2*N )];
   sXF(iiT,:)=imag(fft([YX(iiT,1:N+1) 0  -YX(iiT,N+1:-1:2)]))/sqrt(2*(N+1));
   sVF(iiT,:)=imag(fft([YV(iiT,1:N+1) 0  -YV(iiT,N+1:-1:2)]))/sqrt(2*(N+1));
   Energ(iiT,1:N)=(omegak2(1:N).*(sXF(iiT,2:N+1).^2)+sVF(iiT,2:N+1).^2)/2;
   for J=2:N-1,   % Space loop
     DifY(iiT,J)=Y(iiT,J+1)-Y(iiT,J);
   end
 end

 figure(3);clf
 plot(TiiME,Energ(:,1),TiiME,Energ(:,2),TiiME,Energ(:,3),TiiME,Energ(:,4));
 surf(DifY);  % Space derivative field to show the soliton dynamics 
end %if(myenergy)