Cs=0.5; lambda=1/18; tau=10; n=50;	% set constants, and do 50 steps
t=1:n; c=zeros(1,n); h=zeros(1,n); dt=1:n-1;	% initialize the data arrays
						% this is the loop in which we calculate the time steps
for i=2:n
   c(i)=c(i-1)+(Cs-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 of concentrations
subplot(211)
p=plot(t,c,'r',t,h,'g'); set(p,'LineWidth',2);
xlabel('Time (in years)'); ylabel('Concentration')
set(gca,'FontSize',16,'LineWidth',2);
subplot(212)
p=semilogy(dt,diff(c),'r',dt,diff(h),'g'); set(p,'LineWidth',2);			% and incremental changes with time
xlabel('Time (in years)'); ylabel('Incremental Change')
set(gca,'FontSize',16,'LineWidth',2);

[c(n) h(n)]