% Generates a Finite Differences anisotropic Laplacian
% Dirichlet b.c. on the whole boundary.
% -epsxU_xx - epsyU_yy = f, U=0 on the boundary

function A=cd_anisot(times,epsx,epsy)

disp(['Size of anisotropy in x: ',num2str(epsx)])
disp(['Size of anisotropy in y: ',num2str(epsy)])

init=2; nanr = init^times;  
hstep=1/nanr; h2=hstep*hstep; nall=nanr^2;

%disp('Begin allocating memory...')
lengthA=nanr*nanr*5;
AI  = zeros(lengthA,1);
AJ  = zeros(lengthA,1);
AV  = zeros(lengthA,1);
nextA = 0;
%disp('...end allocating memory.')%

%
AD = spdiags(2*(epsx+epsy)*ones(nall,1),0,nall,nall);

for k=1:nall-1, 
    nextA    = nextA + 1;
    AI(nextA) = k;
    AJ(nextA) = k+1;
    AV(nextA) = -epsx;
end
for k=2:nall,   
    nextA    = nextA + 1;
    AI(nextA) = k;
    AJ(nextA) = k-1;
    AV(nextA) = -epsx;
end
for k=1:nall-nanr, 
    nextA    = nextA + 1;
    AI(nextA) = k;
    AJ(nextA) = k+nanr;
    AV(nextA) = -epsy;
end
for k=nanr+1:nall,   
    nextA    = nextA + 1;
    AI(nextA) = k;
    AJ(nextA) = k-nanr;
    AV(nextA) = -epsy;
end

for k=1:nanr-1,
    ind=k*nanr;
    if ind+1<=nall, 
       nextA    = nextA + 1;
       AI(nextA) = ind;
       AJ(nextA) = ind+1;
       AV(nextA) = epsx;
    end
    if ind+1<=nall, 
       nextA    = nextA + 1;
       AI(nextA) = ind+1;
       AJ(nextA) = ind;
       AV(nextA) = epsx;
    end
end

AI = AI(1:nextA);
AJ = AJ(1:nextA);
AV = AV(1:nextA);
A  = AD + sparse(AI,AJ,AV);

A  = A/(hstep^2);

return
%