function autocatalytic

y0 = [0; 0];
tspan = [0,1000];
opts = odeset('RelTol',1e-4,'AbsTol',1e-10);

delta = 25;
average = input('Average 1 time or 2 times? ');
avg_opt = [average,1];
sol = odeavg(avg_opt,delta,@ode,tspan,y0,opts);

solode = ode23(@ode,tspan,y0,opts);

plot(sol.x,sol.y(1,:),'r',solode.x,solode.y(1,:),'k')
legend('Averaged solution','Computed solution')
title(['Using \Delta = ',num2str(delta),' and average = ',...
        int2str(sol.average),'.'])
xlabel('x(t), which is component 1.')    

figure
plot(sol.x,sol.y(2,:),'r',solode.x,solode.y(2,:),'k')
legend('Averaged solution','Computed solution')
title(['Using \Delta = ',num2str(delta),' and average = ',...
        int2str(sol.average),'.'])
xlabel('y(t), which is component 2.')    

%===Nested function================================================
function yp = ode(t,y)
    yp = [0; 0];
    yp(1) = exp(-0.002*t) - 0.08*y(1) - y(1)*y(2)^2;
    yp(2) = 0.08*y(1) - y(2) + y(1)*y(2)^2;
end % ode
%==================================================================

end % autocatalytic