clear all;


global e f;

% e = van der pol parameter
% f = forcing parameter
e=0.1;
f=0;

% initial conditions
x0=2;
y0=0;

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


% if you set this to =1 then the 
% program will also plot the solution
% from the normal perturbation method (which resonates)
plot_pert=0;

% 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 circle with radius the initial conditions
th = 0:pi/50:2*pi;
xunit = x0 * cos(th);
yunit = x0 * sin(th);



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

Xperturb   = x0.* cos(T) + (3/8*x0^2*e).*(-T.*sin(T));

Xtwotimes  = x0.* cos((1+3/8*x0^2*e)*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] = DUF(x,y,t);
[k2,m2] = DUF(x+k1*dt/2,y+m1*dt/2,t+dt/2);
[k3,m3] = DUF(x+k2*dt/2,y+m2*dt/2,t+dt/2);
[k4,m4] = DUF(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,10)==1
    % plot X(t)
    figure(1),clf;
    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);
    if plot_pert==1
        plot(T(1:it),Xperturb(1:it),'--k','LineWidth',2);
    end
    hold off;
    xlabel('$t$','FontSize',26,'Interpreter','latex');
    ylabel('$x$','FontSize',26,'Interpreter','latex');
    if plot_pert==1
        ylim([-10 10])
    else
        ylim([-5 5])
    end
    if T(it)<80
        xlim([0 100])
        if plot_pert==1
            xlim([0 Tf])
        end
    else
        xlim([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',26,'Interpreter','latex');
    ylabel('$dx/dt$','FontSize',26,'Interpreter','latex');
    ylim([-5 5])
    if T(it)<80
        xlim([0 100])
        if plot_pert==1
            xlim([0 Tf])
        end
    else
        xlim([round(Tf/4) Tf])
    end
    % 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',26,'Interpreter','latex');
    ylabel('$dx/dt$','FontSize',26,'Interpreter','latex');
    r=max([max(X),max(Y),x0,y0,4]);
    axis([-r r -r r]);
    axis square
    drawnow;
    end


end