function x0new = Advection(InVec,CVec,p,dt,nsteps)

% setup the mesh on which to solve
% p=40; m=p^2;
m=p^2;
dx = 1/(p+1);
x0 = InVec;

% ---------------------------------------------------------------------------------------------------------------

% compute the initial solution profile x0 (soln to heat equation)
Ab=zeros(m,p+1); Ab(:,1)=4; Ab(:,2)=-1; Ab(p:p:m,2)=0; Ab(:,p+1)=-1;
A=spdiags(Ab(:,[1,2,p+1]),-[0,1,p],m,m); A=A+tril(A,-1)';
% f=ones(m,1)/(p+1)^2;
% x0=A\f;

% ---------------------------------------------------------------------------------------------------------------

% Calculate the first derivative matrices
Ac_xdirn=zeros(m,p+1); Ac_xdirn(:,1)=1; Ac_xdirn(p+1:p:m+1,1)=0; Ac_xdirn(:,2)=-1; Ac_xdirn(p:p:m,2)=0; Ac_xdirn(:,3)=0; Ac_xdirn(1:p:m-p+1,3)=1; Ac_xdirn(:,4)=0; Ac_xdirn(1:p:m-p+1,4)=-1;
Ac_ydirn=zeros(m,2); Ac_ydirn(:,1)=-1; Ac_ydirn(:,2)=1;
A_yfirstderiv=spdiags(Ac_ydirn(:,[1,1,2,2]),[m-p,-p,p-m,p],m,m);
A_xfirstderiv=spdiags(Ac_xdirn(:,[2,3]),[-1,-p+1],m,m);
A_xfirstderiv = A_xfirstderiv - tril(A_xfirstderiv,-1)';

% disp('full(A_xfirstderiv)')
% full(A_xfirstderiv)
% disp('full(A_yfirstderiv)')
% full(A_yfirstderiv)
% ---------------------------------------------------------------------------------------------------------------

% Calculate the second derivative matrices
d2_xdirn=zeros(m,p+1); d2_xdirn(:,1)=-2; d2_xdirn(p+1:p:m+1,1)=0; d2_xdirn(:,2)=1; d2_xdirn(p:p:m,2)=0; d2_xdirn(:,3)=0; d2_xdirn(1:p:m-p+1,3)=1; d2_xdirn(:,4)=0; d2_xdirn(1:p:m-p+1,4)=-1;
d2_xderiv=spdiags(d2_xdirn(:,[1,2,3]),[0,-1,-p+1],m,m);
d2_xderiv = d2_xderiv - tril(d2_xderiv,-1)';
d2_ydirn=zeros(m,3); d2_ydirn(:,1)=1; d2_ydirn(:,2)=1; d2_ydirn(:,3)=-2;
d2_yderiv=spdiags(d2_ydirn(:,[1,1,2,2,3]),[m-p,-p,p-m,p,0],m,m);

% disp('full(d2_xderiv)')
% full(d2_xderiv)
% disp('full(d2_yderiv)')
% full(d2_yderiv)
% ---------------------------------------------------------------------------------------------------------------

% Compute gradients of 'c' to obtain x and y advection velocities
% do this tomorrow....

ConvectionVel_xdir = A_xfirstderiv*CVec;
% norm(ConvectionVel_xdir)
ConvectionVel_ydir = A_yfirstderiv*CVec;
% norm(ConvectionVel_ydir)

ConvectionVel_xdir = 15*ConvectionVel_xdir;
ConvectionVel_ydir = 15*ConvectionVel_ydir;

% ConvectionVel_xdir = ones(m,1);
% ConvectionVel_ydir = ones(m,1);
% ConvectionVel_xdir
% ConvectionVel_ydir

% alph = 0.01;
% xVel = 0.01;
% yVel = 0.0;
% yVel = 0.01;
% DiffCoeff1x = (xVel.*dt)/(2*dx);
% DiffCoeff1y = (yVel.*dt)/(2*dx);

% alpha = 0.01;
alpha = 1;
DiffCoeff1x = (alpha*dt)/(2*dx);
DiffCoeff1y = (alpha*dt)/(2*dx);

DiffCoeff2x = 2*(DiffCoeff1x.^2);
DiffCoeff2y = 2*(DiffCoeff1y.^2);

% DiffCoeff1 = Vel*dt/(2.*dx);
% DiffCoeff2 = Vel*2*(DiffCoeff1.^2);

% full(A_xfirstderiv)
% full(A_yfirstderiv)
% full(A_xsecondderiv)
% full(A_ysecondderiv)

% I = sparse(1:m,1:m,1);

% ConvectionVel_xdir = 10;
% ConvectionVel_ydir = 10;

for i = 1:nsteps
%     disp('did an advection!')
    x0new = x0 - DiffCoeff1x*A_xfirstderiv*(ConvectionVel_xdir.*x0) - DiffCoeff1y*A_yfirstderiv*(ConvectionVel_ydir.*x0) - DiffCoeff2x*d2_xderiv*(ConvectionVel_xdir.*ConvectionVel_xdir.*x0) - DiffCoeff2y*d2_yderiv*(ConvectionVel_ydir.*ConvectionVel_ydir.*x0);
%     x0new = x0 - DiffCoeff1x*A_xfirstderiv*(ConvectionVel_xdir.*x0) - DiffCoeff1y*A_yfirstderiv*(ConvectionVel_ydir.*x0) - DiffCoeff2x*d2_xderiv*(ConvectionVel_xdir.*ConvectionVel_xdir.*x0) - DiffCoeff2y*d2_yderiv*(ConvectionVel_ydir.*ConvectionVel_ydir.*x0);

%     x0new = (I - DiffCoeff1x*A_xfirstderiv - DiffCoeff1y*A_yfirstderiv - DiffCoeff2x*d2_xderiv - DiffCoeff2y*d2_yderiv)*x0;
%     close;
%     if mod(i,100) == 0
%         PlotData(p,x0new)
%         pause(0.1)
%     end
    x0 = x0new;
end