function fig613

% clear all previous variables and plots
clear *
clf

tmax=1;
ep=0.001; 

% get(gcf)
%set(gcf,'Position', [805 516 578 232]);
set(gcf,'Position', [1896 1238 573 199]);

% steady state
ns=100;
ys=linspace(-3,2,ns);
for i=1:ns
	ts(i)=ys(i)-ys(i)^3/3;
end;

% asy sol
na=200;
taM=2/3+2.338*ep^(2/3);
tam=0.65;
ta=linspace(tam,taM,na);
for i=1:na
	tt=(ta(i)-2/3)/ep^(2/3);
	ya(i)=1-ep^(1/3)*airy(1,-tt)/airy(-tt);
end;

%  initial values
y10=2;  
y0=[y10];

%  calculate solution using a MATLAB routine
[t,y] = ode45(@rhs,[0 tmax],y0);  

hold on

plot(t,y,'--','LineWidth',1.1)
plot(ta,ya,'-r','LineWidth',1)
%plot(ts,ys,':','LineWidth',2)

xp1=2/3; yp1=1;
plot(xp1,yp1,'.k','LineWidth',1,'MarkerSize',24)

axis([0.63 0.73 -3 2]);
% commands to label each axes
xlabel('t-axis','FontSize',12,'FontWeight','bold')
ylabel('Solution','FontSize',12,'FontWeight','bold')
grid on

% command to put legend into plot
loc='SouthWest';
legend(' Numerical',' Corner Layer',' Bifurcation Point','Location',loc);

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

hold off

%  define f(t,y)
function dy=rhs(t,y)
dy=zeros(1,1);
ep=0.001; 
dy(1) = ( y(1) - y(1)^3/3 - t )/ep;