close all
clear all

Nstep = 16;

gamma = 0.5;
omega = 5.0;

u_0 = [0; 1];

t_0 =0.0;
T   = 10.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);

ndiv = 1:Nstep;
Dt = 0.02./(2.^ndiv);
for step = ndiv
  tk = t_0:Dt(step):(t_0+T);
  K = length(tk);
  uk = zeros(2,K);
  uk(:,1) = u_0;

  x  = c1*exp(l1*tk)+c2*exp(l2*tk);

  for k=1:K-1
    uk(:,k+1) = m_theta(0.0,omega,gamma,Dt(step),uk(:,k));
  end
  e_inf_theta_0(step) = norm(abs(uk(1,:)-x),inf);

  for k=1:K-1
    uk(:,k+1) = m_theta(0.5,omega,gamma,Dt(step),uk(:,k));
  end
  e_inf_theta_05(step) = norm(abs(uk(1,:)-x),inf);

  for k=1:K-1
    uk(:,k+1) = m_theta(1.0,omega,gamma,Dt(step),uk(:,k));
  end
  e_inf_theta_1(step) = norm(abs(uk(1,:)-x),inf);

  for k=1:K-1
    uk(:,k+1) = rk4(omega,gamma,Dt(step),uk(:,k));
  end
  e_inf_rk4(step) = norm(abs(uk(1,:)-x),inf);
end


figure(1)
handle = axes
set(handle,'FontSize',18)
handle = loglog(Dt,e_inf_theta_0,'--');
set(handle,'LineWidth',4);
axis([1.0E-7 1.0E-2 1.0E-20 1])
grid on
hold on

handle = loglog(Dt,e_inf_theta_05,'-.');
set(handle,'LineWidth',4);

handle = loglog(Dt,e_inf_theta_1,':');
set(handle,'LineWidth',4);

handle = loglog(Dt,e_inf_rk4,'-');
set(handle,'LineWidth',4);

legend('\theta=0.0','\theta=0.5','\theta=1.0','rk4',4)
xlabel('Dt')
ylabel('error')

hold off

print -deps2 errors.eps