% an example of non-hyperbolic fixed point
% linearized system would yield wrong information on the local behavior for
% the damped oscillator system

tinit=0;
tend=100;
x0=0.5;
y0=0.5;


u0=[x0,y0];

% a simple model for a damped oscillator and ...
[t,sol]= ode23s(@damped_rhs,[tinit tend],u0);

%.. a model for a harmonic oscillator ...
[t2,sol2]= ode23s(@harmon_rhs,[tinit tend],u0);

% ...have the same linearization at the equilibrium point (0,0)
[t3,sol3]= ode23s(@linear_rhs,[tinit tend],u0);

% plot solutions - the origin is a center for the linearized system and for
% the harmonic oscillator, but not for the damped oscillator

figure(1)
subplot(3,3,1)
plot(t,sol(:,1),'b-','LineWidth',1.5)
title('x')
ylabel('DAMPED MODEL')

subplot(3,3,2)
plot(t,sol(:,2),'k-','LineWidth',1.5)
title('y')

subplot(3,3,3)
plot(sol(:,1),sol(:,2),'r-','LineWidth',1.5)
title('Phase space')
xlabel('x')
ylabel('y')


subplot(3,3,4)
plot(t2,sol2(:,1),'b:','LineWidth',1.5)
title('x')
ylabel('HARMONIC OSCIL')

subplot(3,3,5)
plot(t2,sol2(:,2),'k:','LineWidth',1.5)
title('y')

subplot(3,3,6)
plot(sol2(:,1),sol2(:,2),'r:','LineWidth',1.5)
title('Phase space')
xlabel('x')
ylabel('y')

subplot(3,3,7)
plot(t3,sol3(:,1),'b.-','LineWidth',1.2)
title('x')
ylabel('LINEAR')

subplot(3,3,8)
plot(t3,sol3(:,2),'k.-','LineWidth',1.2)
title('y')

subplot(3,3,9)
plot(sol3(:,1),sol3(:,2),'r.-','LineWidth',1.2)
title('Phase space')
xlabel('x')
ylabel('y')


function dy=damped_rhs(t,y,P) 

% --------------------------------------------------------------------
% equations for a damped oscillator
% --------------------------------------------------------------------

dy(1,1)= y(2); %
dy(2,1)= -y(2)*y(1).^2-sin(y(1)); % 
end



function dy=harmon_rhs(t,y,P) 

% --------------------------------------------------------------------
% equations for an harmonic oscillator
% --------------------------------------------------------------------

dy(1,1)= y(2); % 
dy(2,1)= -sin(y(1)); % 
end



function dy=linear_rhs(t,y,P) 

% --------------------------------------------------------------------
% linearization at equilibrium (0,0) 
% of both harmonic and damped oscillator model
% --------------------------------------------------------------------

dy(1,1)= y(2); % 
dy(2,1)= -y(1); %
end



