clear all

global a e
figure(1);clf
figure(2);clf

a=0.24;e=2.;

Nt=100;
T=1;
tt=linspace(0,1,Nt);
dt=tt(2)-tt(1);
x3=zeros(Nt,1);

X=linspace(0,1,31);
p=0*X;

for ix=1:length(X);
     x3(1)=X(ix);
     tt(1)=0;


for it=2:Nt
    %Runge-Kutta
    k1=fdot1d(x3(it-1),tt(it-1));
    k2=fdot1d(x3(it-1)+k1*dt/2,tt(it-1)+dt/2);
    k3=fdot1d(x3(it-1)+k2*dt/2,tt(it-1)+dt/2);
    k4=fdot1d(x3(it-1)+k3*dt,tt(it-1)+dt);
    
    x3(it)=x3(it-1)+(dt/6)*(k1+2*k2+2*k3+k4);
    tt(it)=tt(it-1)+dt;
    
    %if x3(it)<0;break;end;
    
end;

  p(ix)=x3(it);
  
  
  
  figure(1)
   plot(tt,x3,'k','linewidth',2)
   xlabel('t');ylabel('x');title(['a = ',num2str(a),'e = ',num2str(e)])
   hold on
   drawnow

  
  
end;

 xa1=1/2+sqrt(1/4-a);xa2=1/2-sqrt(1/4-a);
 
   figure(2)
   plot(X,p,'r',X,X,'k',X,0*X','--k',xa1,xa1,'ob',xa2,xa2,'ob','linewidth',2)
   ylim([0 1])
   xlabel('x');ylabel('p(x)');title('Poincare map')
   

  figure(1);hold off