function transupwind(cfun,N,T,ic);
% TRANSUPWIND Solves transport equation (4.4) 
% by the upwind method:
%         u_t + c(x,t) u_x = 0
% on interval 0 < x < 10, 0 < t < T.  
% The coefficient function c(x,t) given by 'cfun'.
% Boundary conditions periodic.
% Uses initial condition specified in ic a row vector of length N.
% And dx=1/N.
% Adaptively determines dt by CFL condition
%     (max |c|) dt/dx = .5 < 1.
% Try "cf = inline('x.^2+1','x','t'); 
%      N=40; x=0:10/N:10; 
%      ic=[zeros(1,N/10) (2-x(N/10+2:2*N/10)) zeros(1,8*N/10+1)];
%      transupwind(cf,N,2,ic);"
% to solve   u_t + (x^2+1) u_x = 0.  Compare to exercise 6 on A#6.
% (Ed Bueler 4/14/2002)

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

% insert initial condition
if not(isequal(size(ic),[1 N])), error('Init cond wrong size.'); end;
w(1,1:N)=ic; w(1,N+1)=w(1,1); % "periodize" for appearance

l=1;
while tt<T
   for k=1:N  %build vector of coeffs:
      coeff(k)=feval(cfun,x(k),tt); end;
   B=max(abs(coeff));
   dt=.5*dx/(max([B 1e-4]));
   tt=tt+dt; l=l+1; R=dt/dx;
   tsave(l)=tt;
   if coeff(1) >= 0.0
      w(l,1)=w(l-1,1)-R*coeff(1)*(w(l-1,1)-w(l-1,N));
   else
      w(l,1)=w(l-1,1)-R*coeff(1)*(w(l-1,2)-w(l-1,1));
   end;
   for k=2:N
      if coeff(k) >= 0.0
         w(l,k)=w(l-1,k)-R*coeff(k)*(w(l-1,k)-w(l-1,k-1));
      else
         w(l,k)=w(l-1,k)-R*coeff(k)*(w(l-1,k+1)-w(l-1,k));
      end;
   end;
   w(l,N+1)=w(l,1); % maintain "periodized" appearance
end;

%display it
mesh(x,tsave,w); axis tight;
xlabel('x axis'); ylabel('t axis'); zlabel('u axis');
title(['Soln to transport eqn (4.4) by upwinding with N=' num2str(N)]);
view(37.5,30);
%view(2)  % useful viewpoint