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.25,omega,gamma,Dt(step),uk(:,k));
  end
  e_inf_theta_025(step) = norm(abs(uk(1,:)-x),inf);

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

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

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


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

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

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

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

%legend('\theta=0.25','\theta=0.50001','\theta=0.51','\theta=0.75',4)
legend('\theta=0.25','\theta=0.5+1.0E-5','\theta=0.75',4)
xlabel('Dt')
ylabel('error')

hold off

print -deps2 errors_051.eps