clear all
close all

gamma   = 0.5;
omega   = 5.0;
u_0     = [1; 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);

t_0     = 0.0;
T       = 10.0;
TSPAN   = [t_0 t_0+T];

OPTIONS = odeset('RelTol',1.0E-5)
OPTIONS = odeset(OPTIONS,'AbsTol',[1.0E-7;1.0E-4])
OPTIONS = odeset(OPTIONS,'Stats','on')

[tk,uk] = ode45(@funz,TSPAN,u_0,OPTIONS,omega,gamma);

K = length(tk);

Dtk = tk(2:K)-tk(1:K-1);

figure(1)
handle = axes
set(handle,'FontSize',18)
handle = plot(Dtk);
set(handle,'LineWidth',4);
grid on
xlabel('k')
ylabel('Dtk')

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

err_u  = abs(uk(:,1)-x);
err_up = abs(uk(:,2)-xp);

figure(2)
handle = axes
set(handle,'FontSize',18)
handle = plot(tk,uk(:,1));
set(handle,'LineWidth',4);
grid on
xlabel('t')
ylabel('x')

figure(3)
handle = axes
set(handle,'FontSize',18)
handle = plot(tk,err_u);
set(handle,'LineWidth',4);
grid on
xlabel('t')
ylabel('erru')

figure(4)
handle = axes
set(handle,'FontSize',18)
handle = plot(tk,err_up);
set(handle,'LineWidth',4);
grid on
xlabel('t')
ylabel('errup')