load natrsf.dat			% load in surface water history

Cs=natrsf(:,2);t=natrsf(:,1);	% time in col. 1, trit in col 2

c=zeros(50,1);h=zeros(50,1);		% initialize data
c(1)=0.321; h(1)=0.179;			% starting values in 1950

tau=10; lambda=1/18; n=length(natrsf);	% other factors

for i=2:n
   c(i)=c(i-1)+(Cs(i)-c(i-1))/tau - lambda*c(i-1);
   h(i)=h(i-1)+(0 -h(i-1))/tau + lambda*c(i-1);
 end

			% plot up time history
subplot(211)
p=plot(t,c,'r',t,h,'g'); set(p,'LineWidth',2);
set(gca,'LineWidth',2,'FontSize',16)
xlabel('Year');ylabel('Concentration');
axis([1950 2000 0 8])