function fig621

%  solve IVP using MATLAB routines
%         y' = f(t,y)  with   y(0) = y0        
%  where y = (y1, y2, y3 , ..., yn) is an n-vector


% clear all previous variables and plots
clear *
clf

% get(gcf)
set(gcf,'Position', [1677 1122 665 250]);
hold on

tmax=500;
y0=[0.1 0 2 2];

%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],y0);  
%[t,y] = ode23s(@rhs,[0 tmax],y0); 
z=y(:,4)-y(:,2);

[tt,yy] = ode45(@rhs2,[0 tmax],y0);  
%[t,y] = ode23s(@rhs,[0 tmax],y0); 
zz=yy(:,4)-yy(:,2);

plot(t,z,'k','LineWidth',1.1)
plot(tt,zz,'--k','LineWidth',1.1)


axis([0 500 -0.2 4]);
% commands to label each axes
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('\phi - axis','FontSize',14,'FontWeight','bold')

set(gca,'xtick',[0 100 200 300 400 500]);

grid on
% command to put legend into plot
legend(' a = 1',' a = -1','Location','East');

% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 
% Set legend font to 14/bold                            		
set(findobj(gcf,'tag','legend'),'FontSize',14,'FontWeight','bold'); 

hold off

figure
set(gcf,'Position', [993 1123 666 247]);

hold on
plot(t,y(:,1),'k')
plot(t,y(:,3),'--k')

%plot(tt,yy(:,1),'k')
%plot(tt,yy(:,3),'--k')

axis([0 10 0 2]);
% commands to label each axes
xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('r - axis','FontSize',14,'FontWeight','bold')

%set(gca,'xtick',[0 100 200 300 400 500]);

grid on
% command to put legend into plot
legend(' r_1',' r_2','Location','NorthEast');

% have MATLAB use certain plot options (all are optional)
box on
% Set the fontsize to 14 for the plot
set(gca,'FontSize',14); 




%  define f(t,y)
function dy=rhs(t,y)
global a0
% y(1)=r(1) y(2)=theta(1)  y(3)=r(2)  y(4)=theta(2)
dy=zeros(4,1);
ep=0.01;
a=1; 
b=1;
a1=1; b1=1;  k1=1;
a2=a1;  b2=b1;  k2=k1;

z1=y(3)*cos(y(4)) - y(1)*cos(y(2));  
z2=y(3)*sin(y(4)) - y(1)*sin(y(2));

f1 = -a*z1 - b*z2;  g1 = b*z1-a*z2;
f2 = a*z1 + b*z2; g2=-b*z1+a*z2;

F1=f1*cos(y(2))+g1*sin(y(2));
G1=(g1*cos(y(2))-f1*sin(y(2)))/y(1);

F2=f2*cos(y(4))+g2*sin(y(4));
G2=(g2*cos(y(4))-f2*sin(y(4)))/y(3);

dy(1)=-b1*y(1)*(y(1)^2-k1^2)+ep*F1;
dy(2)=-a1+ep*G1;
dy(3)=-b2*y(3)*(y(3)^2-k2^2)+ep*F2;
dy(4)=-a2+ep*G2;



%  define f(t,y)
function dy=rhs2(t,y)
global a0
% y(1)=r(1) y(2)=theta(1)  y(3)=r(2)  y(4)=theta(2)
dy=zeros(4,1);
ep=0.01;
a=-1; 
b=1;
a1=1; b1=1;  k1=1;
a2=a1;  b2=b1;  k2=k1;

z1=y(3)*cos(y(4)) - y(1)*cos(y(2));  
z2=y(3)*sin(y(4)) - y(1)*sin(y(2));

f1 = -a*z1 - b*z2;  g1 = b*z1-a*z2;
f2 = a*z1 + b*z2; g2=-b*z1+a*z2;

F1=f1*cos(y(2))+g1*sin(y(2));
G1=(g1*cos(y(2))-f1*sin(y(2)))/y(1);

F2=f2*cos(y(4))+g2*sin(y(4));
G2=(g2*cos(y(4))-f2*sin(y(4)))/y(3);

dy(1)=-b1*y(1)*(y(1)^2-k1^2)+ep*F1;
dy(2)=-a1+ep*G1;
dy(3)=-b2*y(3)*(y(3)^2-k2^2)+ep*F2;
dy(4)=-a2+ep*G2;