

clear all

global omega gamma f

omega=1;
f=1;
gamma=1;

N=2;

x0=[1;0];

dt=0.1;

Nt=1000;

xt=x0;
x1(1)=xt(1);
x2(1)=xt(2);

t(1)=0;

for it=2:Nt;

    t(it)=(it-1)*dt;
    tt=t(it);
    %Euler
    %k1=dt*vel2d(xt,tt);
    %xt=xt+k1;
    %RK4
    
    k1=dt*vel2d(xt,tt);
    k2=dt*vel2d(xt+k1/2,tt+dt/2);
    k3=dt*vel2d(xt+k2/2,tt+dt/2);
    k4=dt*vel2d(xt+k3,tt+dt);
    xt=xt+(k1+2*k2+2*k3+k4)/6;


    x1(it)=xt(1);
    x2(it)=xt(2);

end
%xan=cos(t);

figure(103)
plot(x1,x2,'r','Linewidth',2)
%plot(t,x,'','Linewidth',2)
xlabel('$x_1$','Interpreter','latex')
ylabel('$x_2$','Interpreter','latex')
set(gca,'Fontsize',24)
grid on
%hold on
figure(104)
plot(t,x1,'r',t,0*x1,'k','Linewidth',2)
%plot(t,x1,'r',t,0*x1,'k',t,xan,'.b','Linewidth',2)
%plot(t,x,'','Linewidth',2)
xlabel('$t$','Interpreter','latex')
ylabel('$x(t)$','Interpreter','latex')
set(gca,'Fontsize',24)
grid on
%hold on


