% Solve the homework problem using a grid of n points
% starting from initial condition IC (function handle)
function [w,x,h] = AdvDiff(a, D, L, t, nu, n, IC, DBC)
   h=L/n;
   
   global periodic

   if a<=0
      dt=h; % Diffusion only
   else  
      dt=nu*h/a; % Time step size
   end   
   n_steps=round(t/dt)
   dt=t/n_steps

   x = h*[1/2:(n-1/2)]';
   x_mid = x+h/2;
   d_mid = D(x_mid); % Diffusion coefficient evalued at faces

   mu=max(d_mid)*dt/h^2 % Diffusive CFL number

   % Form the diffusive operator:
   %-----------------------
   if(periodic) 
      d_mid_m = circshift(d_mid, 1);
      A = spdiags([d_mid /h^2 , -(d_mid + d_mid_m)/h^2, circshift(d_mid,1)/h^2],[-1 0 1],n,n);
      % Correct corners for periodic boundary conditions:
      A(1, n) = d_mid_m(1)/h^2;
      A(n,1) = d_mid(n)/h^2;
   else
      d_mid_m = circshift(d_mid, 1);
      A = spdiags([d_mid /h^2 , -(d_mid + d_mid_m)/h^2, circshift(d_mid,1)/h^2],[-1 0 1],n,n);
      % Correct first and last row for Dirichlet/Neumann BC
      A(1,1) = -(2*d_mid(1)+d_mid(2))/h^2; A(1,2) = d_mid(2)/h^2;
      A(n,n-1) = d_mid(n-1)/h^2; A(n,n) = -d_mid(n-1)/h^2; 
   end

   % Time stepping:
   %-----------------------
   w0=IC(x); % Initial condition
   w=w0;
   for step=1:n_steps
      w = TimeStep(w, A, a, dt, h, n, (step-1)*dt, d_mid, DBC); 
         % Matlab allows you to have w be both input and output
         % But in other programming languages it may be different
   end

end
