function fig19

clf

nt=200;
t=linspace(0,2,nt);

% get(gcf)
set(gcf,'Position', [1041 771 548 229]);

ep=0.1;
for it=1:nt
	ya1(it)=t(it)*(1-0.5*t(it));
	ya2(it)=ya1(it)+ep*t(it)^3*(1-0.25*t(it))/3;
end;

%  exact
y10=0; y20=1; 
y0=[y10 y20];
[tt,yy] = ode45(@rhs,[0 2],y0);  
%[t,y] = ode23s(@rhs,[0 2],y0); 

hold on

plot(tt,yy(:,1),'-','Linewidth',1)
plot(t,ya1,'--','Linewidth',1)
plot(t,ya2,'-.','Linewidth',1)

box on
grid on
axis([0 2 0 0.6])
%loc='NorthWest';
loc='South';

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

set(gca,'ytick',[0 0.2 0.4 0.6]);
set(gca,'FontSize',14);
legend(' Numerical Solution',' y_0',' y_0+\epsilon y_1','Location',loc);
set(findobj(gcf,'tag','legend'),'FontSize',14); 


%  define f(t,y)
function dy=rhs(t,y)
dy=zeros(2,1);
ep=0.1; 
dy(1) = y(2);
dy(2) = -1/(1+ep*y(1))^2;