function [x,t,N]=AdvDiff1D(N0,u,k,L,T,dh,dt);
%ADVDIFF1D--simulates simple 1-D advection-diffusion equations
%
%[x,t,N]=AdvDiff1D(N0,u,k,L,T,dh,dt);
%
% N0=initial density.  Initial dist. is N0*(sin(pi*x/L).^2);
% u=velocity
% k=diffusion coefficient
% L=Length of domain
% T=max time of simulation
% dh=grid spacing
% dt=time step

tmp=L/dh;
m=round(tmp);
dh=L/m;
x=[0:dh:L]';m=m-1;
% remove endpoints for periodic BC's
tmp=T/dt;
n=round(tmp);
dt=T/n;
t=[0:dt:T];

sigma=k*dt/dh^2;
lambda=u*dt/dh;
disp(['sigma=',num2str(sigma),', lambda=',num2str(abs(lambda))]);
nnz=3*(m-2)+4; %number of non-zero entries in the matrix
I=zeros(nnz,1);
S=I;
%main diagonal
I(1:m)=[1:m]';
J=I;
S(1:m)=1+sigma;
%upper diagonal
I(m+1:2*m-1)=[1:m-1]';
J(m+1:2*m-1)=[2:m]';
S(m+1:2*m-1)=-0.5*sigma;
%lower diagonal
I(2*m:3*m-2)=[2:m]';
J(2*m:3*m-2)=[1:m-1]';
S(2*m:3*m-2)=-0.5*sigma;
%load Matrix;
M=sparse(I,J,S);
M(1,m)=-0.5*sigma;%periodic BC's
M(m,1)=-0.5*sigma;
%factor
[LT,UT]=lu(M);

N=zeros(m+2,n);

N(:,1)=N0*(sin(pi*x/L).^2);
%I=round(0.25*(m+2)):round(0.75*(m+2));
%N(I,1)=N0;


for j=2:n;
	%RHS=N(2:m+1,j-1)*(1+dt*r)+1/2*sigma*(N(3:m+2, j-1)-2*N(2:m+1,j-1)...
  	%																			+N(1:m,j-1));
   N(1,j-1)=N(m+1,j-1);%periodic BCs
   N(m+2,j-1)=N(2,j-1);
   RHS=N(2:m+1,j-1)+0.5*sigma*( N(3:m+2, j-1)-2*N(2:m+1,j-1)...
  																				+N(1:m,j-1) );
   %add Lax-Wendroff to RHS
   RHS=RHS+ 0.5*lambda*(N(3:m+2, j-1)-N(1:m,j-1));
   RHS=RHS+ 0.5*lambda*lambda*( N(3:m+2, j-1)-2*N(2:m+1,j-1) +N(1:m,j-1) );
   			
  RHS=LT\RHS;
  N(2:m+1,j)=UT\RHS;
end
N([1,m+2],:)=[];
x([1,m+2])=[];