function fig33

clf

tmax=150;

% get(gcf)
set(gcf,'Position', [1665 1089 638 257]);

ep=0.05; lambda=2; w=3/8; kap=4;

%  exact
y10=0; y20=0; 
y0=[y10 y20];
[t,y] = ode45(@rhs1,[0 tmax],y0);  
%[t,y] = ode23s(@rhs1,[0 tmax],y0); 

%  asy
y10=0.000000001; y20=-0.5*pi; 
y0=[y10 y20];
etmax=ep*tmax;
[tt,yy] = ode45(@rhs2,[0 etmax],y0);  
%[tt,yy] = ode23s(@rhs2,[0 etmax],y0);
tt=tt/ep; 

hold on
box on
grid on

plot(t,y(:,1),'-r','Linewidth',1)
plot(tt,yy(:,1),'--','Linewidth',1.3)

axis([0 tmax -0.601 0.6])
%loc='NorthWest';
loc='SouthEast';

xlabel('t-axis','FontSize',14,'FontWeight','bold')
ylabel('Solution','FontSize',14,'FontWeight','bold')

set(gca,'FontSize',14);
%legend(' Exact Solution',' y_0(t)+\epsilon y_1(t)','Location',loc);
set(findobj(gcf,'tag','legend'),'FontSize',14); 


%  define f(t,y)
function dy=rhs1(t,y)
dy=zeros(2,1);
ep=0.05; lambda=2; w=3/8; kap=4;
ff=ep*cos((1+ep*w)*t);
dy(1) = y(2);
dy(2) = -y(1)-ep*lambda*y(2)-ep*kap*y(1)^3+ff;


%  define f(t,y)
function dy=rhs2(t,y)
dy=zeros(2,1);
ep=0.05; lambda=2; w=3/8; kap=4;
dy(1) = 0.5*( -lambda*y(1) - sin(y(2)-w*t));
dy(2) = 0.5*(0.75*kap*y(1)^3-cos(y(2)-w*t))/y(1);