close all
clear all

gamma = 0.5;
omega = 5.0;

u_0 = [0; 1];

Dt = 0.01;

t_0 =0.0;
T   = 10.0;

tk = t_0:Dt:(t_0+T);

K = length(tk);

uk = zeros(2,K);

uk(:,1) = u_0;

l1 = -gamma+sqrt(gamma^2-omega^2);
l2 = -gamma-sqrt(gamma^2-omega^2);
c1 = (u_0(2)-l2*u_0(1))/(l1-l2);
c2 = -(u_0(2)-l1*u_0(1))/(l1-l2);
x  = c1*exp(l1*tk)+c2*exp(l2*tk);
xp = l1*c1*exp(l1*tk)+l2*c2*exp(l2*tk);

figure(1)
handle = axes
set(handle,'FontSize',18)
handle = plot(tk,x,'g:');
set(handle,'LineWidth',4);
axis([0 10 -0.14 0.20])
grid on

hold on

for k=1:K-1
  uk(:,k+1) = m_theta(0.5,omega,gamma,Dt,uk(:,k));
end
handle = plot(tk,uk(1,:),'b-');
set(handle,'LineWidth',2);

for k=1:K-1
  uk(:,k+1) = m_theta(0.75,omega,gamma,Dt,uk(:,k));
end
handle = plot(tk,uk(1,:),'b-.');
set(handle,'LineWidth',2);

for k=1:K-1
  uk(:,k+1) = m_theta(1.0,omega,gamma,Dt,uk(:,k));
end
handle = plot(tk,uk(1,:),'b--');
set(handle,'LineWidth',2);

legend('esatta','\theta=0.5','\theta=0.75','\theta=1.0')

hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(2)
handle = axes
set(handle,'FontSize',18)
handle = plot(tk,x,'g:');
set(handle,'LineWidth',4);
%axis([0 10 -4 4]) %Dt = 0.1
axis([0 10 -0.14 0.20]) %Dt = 0.01
grid on
hold on

for k=1:K-1
  uk(:,k+1) = m_theta(0.0,omega,gamma,Dt,uk(:,k));
end
handle = plot(tk,uk(1,:),'b-.');
set(handle,'LineWidth',2);

for k=1:K-1
  uk(:,k+1) = m_theta(0.25,omega,gamma,Dt,uk(:,k));
end
handle = plot(tk,uk(1,:),'b--');
set(handle,'LineWidth',2);

legend('esatta','\theta=0.0','\theta=0.25',3)

hold off

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(3)
handle = axes
set(handle,'FontSize',18)
handle = plot(x,xp,'g:');
set(handle,'LineWidth',4);
axis([-0.2 0.2 -1 1])
grid on
hold on

for k=1:K-1
  uk(:,k+1) = m_theta(0.0,omega,gamma,Dt,uk(:,k));
end
handle = plot(uk(1,:),uk(2,:),'b-.');
set(handle,'LineWidth',2);

for k=1:K-1
  uk(:,k+1) = m_theta(0.5,omega,gamma,Dt,uk(:,k));
end
handle = plot(uk(1,:),uk(2,:),'b');
set(handle,'LineWidth',2);

for k=1:K-1
  uk(:,k+1) = m_theta(1.0,omega,gamma,Dt,uk(:,k));
end
handle = plot(uk(1,:),uk(2,:),'b--');
set(handle,'LineWidth',2);

legend('esatta','\theta=0.0','\theta=0.5','\theta=1.0',3)

xlabel('x')
ylabel('xp')

hold off