function A = advecdiff2(N,a,d)
% Central finite difference discretization of
% - div (D grad u)  + < a,grad u > = f
% with Dirichlet zero boundary conditions (?)
% Input:
%     N: Number of interior grid points (mesh spacing 1/(N+1))
%     a: function of x and y with two components 
%     d: function of x and y with two components which are the diagonal
%        of the diffusion coefficient matrix D = [d_1, 0; 0 d_2] 

h = 1/(N+1);
x = linspace(0,1,N+2); x = x(2:end-1);
y = x;
[X,Y] = meshgrid(x,y); 
C = a(X,Y);
C1 = C(1:end/2,:);
C1 = reshape(C1,numel(C1),1);
C2 = C(end/2+1:end,:);
C2 = reshape(C2,numel(C2),1);

D = d(X,Y);
D1 = D(1:end/2,:);
D1 = reshape(D1,numel(D1),1);
D2 = D(end/2+1:end,:);
D2 = reshape(D2,numel(D2),1);


e = ones(N,1);

%A1Grad_neg  = spdiags([-e e],[-1,0],N,N);
%A1Grad_pos  = spdiags([-e e],[ 0,1],N,N);
A1Delta = spdiags([-e 2*e -e],[-1 0 1],N,N);
A1Advec = spdiags([-e e],[-1 1],N,N);

N2 = N*N;

A2Delta = spdiags(D2,0,N2,N2)*kron(speye(N),A1Delta) + spdiags(D1,0,N2,N2)*kron(A1Delta,speye(N));
A2Advec = spdiags(C2,0,N2,N2)*kron(speye(N),A1Advec) + spdiags(C1,0,N2,N2)*kron(A1Advec,speye(N));

A = (1/h)^2*A2Delta + (1/(2*h))*A2Advec;




