clear all;


global e f;

% e = van der pol parameter
% f = forcing parameter
e=.1;
f=0;

% initial conditions
x0=4;
y0=0;

% time step
dt=.1;
Tf=200;

% time vector
T = [0:dt:Tf];

X=0*T;
Y=0*T;
% inserting initial conditions in beginning
% of X and Y vectors
X(1)=x0;
Y(1)=y0;


% plot the limit cycle
th = 0:pi/50:2*pi;
xunit = 2 * cos(th);
yunit = 2 * sin(th);



r0 = sqrt(x0^2+y0^2);

%the analytic solution from the two-time perturbation method
Xtwotimes = ( 2./(sqrt(1 - (1-4/r0^2)*exp(-e*T)) )).* cos(T);



figure(1),clf;



% do the Runge-Kutta integration
for it=2:length(T);

    x=X(it-1);
    y=Y(it-1);
    t=T(it-1);
    [k1,m1] = VP(x,y,t);
    [k2,m2] = VP(x+k1*dt/2,y+m1*dt/2,t+dt/2);
    [k3,m3] = VP(x+k2*dt/2,y+m2*dt/2,t+dt/2);
    [k4,m4] = VP(x+k3*dt,y+m3*dt,t+dt);

    X(it)=x+(k1+2*k2+2*k3+k4)*dt/6;
    Y(it)=y+(m1+2*m2+2*m3+m4)*dt/6;


    if rem(it,round(2/dt))==1
    % plot X(t)
    figure(1);
    subplot(2,4,[1:2])
    plot(T(1:it),X(1:it),'r','LineWidth',2);
    hold on;
    plot(T(1:it),Xtwotimes(1:it),'--b','LineWidth',2);
    hold off;
    xlabel('$t$','FontSize',18,'Interpreter','latex');
    ylabel('$x$','FontSize',18,'Interpreter','latex');
    title('red: RK4, blue: two-time analytic solution','FontSize',18,'Interpreter','latex');
    ylim([-3 3])
    ylim([-10 10])
    if T(it)<80
        xlim([0 100])
    else
        xlim([0*round(Tf/4) Tf])
    end
    

    % plot Y(t)
    subplot(2,4,[5:6])
    plot(T(1:it),Y(1:it),'r','LineWidth',2)
    xlabel('$t$','FontSize',18,'Interpreter','latex');
    ylabel('$dx/dt$','FontSize',18,'Interpreter','latex');
    ylim([-5 5])
    xlim([0 Tf])
    % plot phase space
    subplot(2,4,[3,4,7,8])
    plot(X(1:it),Y(1:it),'r','LineWidth',2);
    hold on;
    plot(xunit,yunit,'--m','LineWidth',2);
    plot(X(it),Y(it),'*r','MarkerSize',10,'LineWidth',2);
    hold off;
    xlabel('$x$','FontSize',18,'Interpreter','latex');
    ylabel('$dx/dt$','FontSize',18,'Interpreter','latex');
    r=1.1*max([max(X),max(Y),x0,y0,4]);
    axis([-r r -r r]);
    axis square
%     drawnow;
    end


end