function burgupwind(c,N,T,ic);
% BURGUPWIND Solves inviscid Burger's equation by the upwind method:
%         u_t = c u u_x
% on interval 0 < x < 1, 0 < t < T.  
% Boundary conditions: u(0,t)=1, u(1,t)=0. (Overdetermines, in fact.)
% Uses initial condition specified in ic a row vector of length N-1.
% And dx=1/N.
% Adaptively determines dt by CFL condition
%     (max |c u|) dt/dx = .5 < 1.
% Try "N=20; x=0:1/N:1; 
%      ic=[ones(1,N/4) (1.5-2*x(N/4+2:3*N/4)) zeros(1,N/4)];
%      burgupwind(-1,N,2,ic);"
% and compare
%      "burgmol(-1,0,N,2,ic,1);"
% (Ed Bueler 4/14/2002)

dx=1/N; x=0:dx:1;
w=zeros(1,N+1); % will grow as needed
tt=0.0; blowup=10.0;
tsave=[tt]; % will record actual t at each step and grow

% insert initial condition
if not(isequal(size(ic),[1 N-1])), error('Init cond wrong size.'); end;
w(1,2:N)=ic;
w(1,1)=1.0;w(1,N+1)=0.0; % boundary cond

l=1;
while tt<T
   B=max(abs(w(l,:))); % >= 1.0 because of bdry cond
   dt=.5*dx/(B*max([abs(c) 1e-4]));
   tt=tt+dt; l=l+1; R=dt/dx;
   tsave(l)=tt; 
   w(l,1)=1.0; % boundary condition
   for k=2:N
      coeff=c*w(l-1,k);
      if coeff <= 0.0
         w(l,k)=w(l-1,k)+R*coeff*(w(l-1,k)-w(l-1,k-1));
      else
         w(l,k)=w(l-1,k)+R*coeff*(w(l-1,k+1)-w(l-1,k));
      end;
   end;
   w(l,N+1)=0;
   if norm(w(l,:),inf) > blowup % stops if size of solution blows up
      error(['blow up at step ' num2str(k)]), end;
end;

%display it
mesh(x,tsave,w); axis tight;
xlabel('x axis'); ylabel('t axis'); zlabel('u axis');
title(['Soln to inviscid Burgers eqn by upwinding for c=' ...
      num2str(c) ' and N=' num2str(N)]);
view(37.5,30);