clear all;


global e f;

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

% initial conditions
x0=1;
y0=0;

% time step
dt=.1;

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

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

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

x=X(it-1);
y=Y(it-1);
t=T(it-1);
[k1,m1] = FG(x,y,t);
[k2,m2] = FG(x+k1*dt/2,y+m1*dt/2,t+dt/2);
[k3,m3] = FG(x+k2*dt/2,y+m2*dt/2,t+dt/2);
[k4,m4] = FG(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;

end
 

% plot X(t)
figure(1)
clf;
plot(T,X,'r');
xlabel('t');
ylabel('x');

% plot Y(t)
figure(2)
plot(T,Y,'r')
xlabel('t');
ylabel('dx/dt');

% plot phase space
figure(3)
plot(X,Y,'b')
xlabel('x');
ylabel('dx/dt');
r=max([max(X),max(Y),x0,y0,4]);
axis([-r r -r r]);
axis square