clear all;

% ''pagkosmies'' statheres
global gamma om

gamma = 0.1;
om    = 1;


% arxikes synthikes
x0 =1;
y0 =0;


% orismos ton xronikon stigmon
T = linspace(0,30,201);
dt = T(2)-T(1)

% to dianysma sto opoio tha apothikevo ti lysi x(t),y(t)
X = zeros(2,length(T));

% i proti colona tou X einai oi arxikes synthikes
X(1,1) = x0;
X(2,1) = y0;



for it=2:length(T)

    x = X(1,it-1);
    y = X(2,it-1);
    t = T(it-1);

    % RK-4
    [k1x , k1y ] = func2D(x,y,t);
    [k2x , k2y ] = func2D(x+k1x*dt/2 , y+k1y*dt/2 , t+dt/2);
    [k3x , k3y ] = func2D(x+k2x*dt/2 , y+k2y*dt/2 , t+dt/2);
    [k4x , k4y ] = func2D(x+k3x*dt  , y+k3y*dt  , t+dt);

    xnew = x + (k1x+2*k2x+2*k3x+k4x)*dt/6;
    ynew = y + (k1y+2*k2y+2*k3y+k4y)*dt/6;

    % apothikevo tis nees times x,y stin epomeni kolona tou X
    X(1,it)=xnew;
    X(2,it)=ynew;

end

% i analytiki lysi gia ton talantoti me aposvesi
OM = sqrt(om^2-gamma^2);
xtheor = x0*exp(-gamma*T).*( cos(OM*T) + gamma/OM*sin(OM*T) ) + y0*exp(-gamma*T)/OM.*sin(OM*T);



figure(1)
plot(T,X(1,:),'ob','LineWidth',2);
hold on;
plot(T,xtheor,'--r')
hold off;
xlabel('t','fontsize',20);
ylabel('x','fontsize',20);
h=legend('numerical','analytical');

figure(3)
plot(X(1,:),X(2,:),'-k')
axis square;
xlabel('x','fontsize',20);
ylabel('xdot','fontsize',20);
drawnow;
pause(1)