 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function exe_cosx
% Nonlinear system:
% 	x1d=x2
% 	x2d=cos(x1)-3*x1/(2*pi)-beta*x2
MainString = mfilename; 
Stampa=0;
global beta
ii=0;
for beta=[-0.5 0:2:4];                      % Beta values
    ii=ii+1;
    figure(ii); clf
    x10=pi/3;                               % Equilibrium point
    x20=0;
    V=[[-1.5 1.5]+x10 [-1.5 1.5]+x20];      % Plot window
    if beta>0
        In_Con=inicond(V,[5,5]);			% Initial conditions
        Tspan=(0:0.01:1)*10;				% Simulation time
        fr=4;								% Arrow's position
    elseif beta==0
        col=x10+(V(2)-x10)*[0.1:1/10:1]';
        In_Con=[col col*0 ];                % Initial conditions
        Tspan=(0:0.01:1)*10;				% Simulation time
        fr=10;                              % Arrow's position
    else
        xx=0.1*exp(1i*(0:0.1:1)*2*pi)';
        In_Con=[real(xx)+x10 imag(xx)+x20]; % Initial conditions
        Tspan=(0:0.01:1)*10;                % Simulation time
        fr=80;                              % Arrow's position
    end
    dx=0.07; dy=dx;
    for jj=(1:size(In_Con,1))
        [t,x]=ode45(@exe_cosx_ode,Tspan,In_Con(jj,:));		% Simulation
        plot(x(:,1),x(:,2)); hold on                        % Data plot
        freccia(x(fr,1),x(fr,2),x(fr+1,1),x(fr+1,2),dx,dy)  % Arrows
    end
    grid on
    axis(V)
    xlabel('Variable x_1(t)')
    ylabel('Variable x_2(t)')
    title(['beta=' num2str(beta)])
    if Stampa;  eval(['print -depsc ' MainString '_' num2str(gcf) '.eps']);  end
end

 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function dx=exe_cosx_ode(t,x)
global beta
dx(1,1)=x(2);
dx(2,1)=cos(x(1))-3*x(1)/(2*pi)-beta*x(2);
return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function  X0=inicond(V, dxy)
% Provides the initial conditions located on a rectangular shape
% V=[x1 x2 y1 y2] with a number of points given by dxy=[dx dy]
 X0=zeros(2*(dxy(1)+dxy(2)),2);
 for ii= 1:dxy(2);  X0(ii,:)=[V(1) V(3)+(ii-1)*(V(4)-V(3))/dxy(2)]; end
 for ii= 1:dxy(1);  X0(dxy(2)+ii,:)=[V(1)+(ii-1)*(V(2)-V(1))/dxy(1) V(4)]; end
 for ii= 1:dxy(2);  X0(dxy(1)+dxy(2)+ii,:)=[V(2) V(4)-(ii-1)*(V(4)-V(3))/dxy(2)]; end
 for ii= 1:dxy(1);  X0(dxy(1)+2*dxy(2)+ii,:)=[V(2)-(ii-1)*(V(2)-V(1))/dxy(1) V(3)]; end
 return

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function f = freccia(x0,y0,x,y,rx,ry)
% Plot an arrow in (x,y) with directio dx=x-x0, dy=y-y0. 
dx=(x-x0)/rx;
dy=(y-y0)/ry;   
dxy=sqrt(dx^2+dy^2);
fx=-dx/dxy;
fy=-dy/dxy;
rotpiu=[rx,0;0,ry]*[cos(pi/6), -sin(pi/6); sin(pi/6), cos(pi/6) ]*[fx,fy]';
rotmeno=[rx,0;0,ry]*[cos(pi/6), sin(pi/6); -sin(pi/6), cos(pi/6) ]*[fx,fy]';
plot([x0,x],[y0,y],'-')
plot([x,x+rotpiu(1)],[y,y+rotpiu(2)],'-')
plot([x,x+rotmeno(1)],[y,y+rotmeno(2)],'-')
return