%Generation of stable limit cycles
%   kja 070820

aminit;

ts=20;tspan=[0,ts];
x01=[0;2];x02=[1 1]/sqrt(2);
[t1,x1]=ode45('limcyc',tspan,x01);
[t2,x2]=ode45('limcyc',tspan,x02);

%   Plot phaseplane
figure(1); clf; subplot(221)
pl=plot(x1(:,1),x1(:,2),'r--');hold on;
set(pl,'Linewidth', AM_data_linewidth);
plot(0, 2, 'ro', 'MarkerSize', AM_data_markersize);
pl=plot(x2(:,1),x2(:,2),'b-');hold on;
set(pl,'Linewidth', AM_data_linewidth);
plot(0.707, 0.707, 'bo', 'MarkerSize', AM_data_markersize);
amaxis([-1.2 2.2 -1.2 2.2]);
line([-1.2 2.2],[0 0], 'LineWidth', AM_ref_linewidth);
line([0 0],[-1.2 2.2], 'LineWidth', AM_ref_linewidth);
axis square; xlabel('x1'); ylabel('x2', 'Rotation', 0);
amprint('limcyc1.eps');

%   Plot time functions
figure(2); clf; subplot(421)
pl=plot(t1,x1(:,1),'r--',t2,x2(:,1),'b-');
set(pl,'Linewidth', AM_data_linewidth);
amaxis([0 20 -1.2 2]);
ylabel('x1', 'Rotation', 0)

subplot(423)
pl=plot(t1,x1(:,2),'r--',t2,x2(:,2),'b-');
set(pl,'Linewidth', AM_data_linewidth);
amaxis([0 20 -1.2 2]);
xlabel('t'); ylabel('x2', 'Rotation', 0);

amprint('limcyc2.eps');