%Perturbation_VanderPol.m

clear all
clc

%options = odeset('OutputFcn','odephas2');

%Tfinal=input('Enter final time:  ');
Tfinal=3;

TSPAN=[0 Tfinal];

%epsilon=input('Enter epsilon value:  ');
epsilon=0.1;

x0=[1;0];

[T,Y]=ode45(@(t,y) Perturbation_VanderPol_fun(t,y,epsilon),TSPAN,x0);

[T0,Y0]=ode45(@(t,y) Perturbation_VanderPol_fun0(t,y,epsilon),T,x0);

[T1,Y1]=ode45(@(t,y) Perturbation_VanderPol_fun1(t,y,epsilon),T,[0;0]);

figure(1)
clf;
plot(T,Y(:,1),'r-',T0,Y0(:,1),'b:',T1,Y0(:,1)+epsilon*Y1(:,1),'k--')
xlabel('Time [s]')
ylabel('x_1')
legend('Exact','0th-order','1st-order','fontsize',14)

figure(2)
clf;
plot(T,Y(:,2),'r-',T0,Y0(:,2),'b:',T1,Y0(:,2)+epsilon*Y1(:,2),'k--')
xlabel('Time [s]')
ylabel('x_2')
legend('Exact','0th-order','1st-order','fontsize',14)


