%		dostep.m
%		does a first upwind differencing method
%		calculation of pure advection of a step
%		function in a linear domain
%
tyr=24*3600*365.25; dx=50000; dt=200000; u=.02;
nt=15*tyr/dt;nx=500;
%			FUDM weights
w0=1-u*dt/dx;
wm=u*dt/dx;
%			Initialize arrays
c=zeros(nx,1); cold=c; x=dx*[1:nx]/1.e6;
oldnow=0;
for i=1:nt
	now=ceil(i*dt/tyr);
	if now ~= oldnow
		subplot(2,2,1)
		p=plot(x,c,'r'); set(p,'LineWidth',2); hold on
		xlabel('Distance');ylabel('Concentration');
		subplot(2,2,2);
		p=plot(x,c,'r'); set(p,'LineWidth',2);
		xlabel('Distance');ylabel('Concentration');		
		pause(1);
		oldnow=now;
	end
c(1)=1;
c(2:nx)=cold(2:nx)*w0 + cold(1:nx-1)*wm;
cold=c;
end
subplot(2,2,1); set(gca,'LineWidth',2,'FontSize',16);
axis([0 12 0 1])
subplot(2,2,2); set(gca,'LineWidth',2,'FontSize',16);
axis([7 11 0 1])