function [t,y] = autocatalyticf90

% Autocatalytic oscillations, Borelli & Coleman p. 260

clc

tempavg = load('autocatalyticf90.dat','-ascii');
tempd = load('autocatalyticf90d.dat','-ascii');
tempode = load('autocatalyticf90ode.dat','-ascii');

delta = tempd(1);
navg = tempd(2);

tode = tempode(:,1);
xode = tempode(:,2);
yode = tempode(:,3);
tavg = tempavg(:,1);
xavg = tempavg(:,2);
yavg = tempavg(:,3);

figure
plot(tavg,xavg,'r',tode,xode,'k')
%legend('Averaged solution','Computed solution')
%title(['Using \Delta = ',num2str(delta),' and average = ',...
%        int2str(navg),'.'])
xlabel('t')    

figure
plot(tavg,xavg,'r',tode,xode,'k')
legend('Averaged solution','Computed solution')
title(['Using \Delta = ',num2str(delta),' and average = ',...
        int2str(navg),'.'])
ylabel('x(t), which is component 1.')    

figure
plot(tavg,yavg,'r',tode,yode,'k')
legend('Averaged solution','Computed solution')
title(['Using \Delta = ',num2str(delta),' and average = ',...
        int2str(navg),'.'])
ylabel('y(t), which is component 2.')