clear all

global a e om

a=1;
e=10;
om=1;


dt=0.01;
t(1)=0;

Tmax=50*2*pi;
N=ceil(Tmax/dt);


Xinit=[0.3;0.1];


    
    x(:,1)=Xinit;
    t(1)=0;
    x1(1)=1;
    x2(1)=0;


for ia=1:N
    
t(ia+1)=ia*dt;
xn=x(:,ia);tn=t(ia);
y = rk_2d(xn,tn,dt);
x(:,ia+1)=y;
x1(ia+1)=0.1*cos(t(ia+1));
x2(ia+1)=-0.1*sin(t(ia+1));
if rem(ia,100)==0
figure(4)
plot(x(1,1:ia+1),x(2,1:ia+1),'k',x1(1:ia+1),x2(1:ia+1),'--r')
xlabel('x_1');ylabel('x_2')
drawnow
figure(5)
plot(t(1:ia+1),x(1,1:ia+1),'b',t(1:ia+1),0*exp(0*t(1:ia+1)),'k')
xlabel('t');ylabel('x_1')
drawnow

end;

end;