clear all;

global a d e omega

%ddot{x}+e*dot(x)*(x^2-1) + x =  a cos (omega t)

e=0.1;
a=1;
omega=0.1;



Tmax=500;


nt=floor(Tmax/0.01);
tt=linspace(0,Tmax,nt);
dt=tt(2)-tt(1);

x00=0.5;y00=0.5;



x0=x00;y0=y00;
%nt=10;


xt(1)=x0;
yt(1)=y0;


for it=2:nt;
    
    tp=tt(it-1);
    %Runge-kutta step
    
    %if tp>100;a=2;end;
    
    xp=xt(it-1);yp=yt(it-1);
    [k1x,k1y]=fvdpol(tp,xp,yp);
    [k2x,k2y]=fvdpol(tp+dt/2,xp+k1x*dt/2,yp+k1y*dt/2);
    [k3x,k3y]=fvdpol(tp+dt/2,xp+k2x*dt/2,yp+k2y*dt/2);
    [k4x,k4y]=fvdpol(tp+dt,xp+k3x*dt,yp+k3y*dt);
    xt(it)=xp + (k1x+2*k2x+2*k3x+k4x)*dt/6;
    yt(it)=yp + (k1y+2*k2y+2*k3y+k4y)*dt/6;
    

    
end;
  
  
figure(1)  
plot(xt,yt,'c')
  
  figure(3)
  plot(tt,xt,tt,0*xt)