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

function a0=matgen_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;

diago = 2*(epsx+epsy);
%
a0=spalloc(nall,nall,5*nall);
for k=1:nall,
    a0(k,k)=diago;
end
for k=1:nall-1, a0(k,k+1)=-epsx;end
for k=2:nall,   a0(k,k-1)=-epsx;end
for k=1:nall-nanr, a0(k,k+nanr)=-epsy; end
for k=nanr+1:nall,   a0(k,k-nanr)=-epsy; end

for k=1:nanr-1,
    ind=k*nanr;
    if ind+1<=nall, a0(ind,ind+1)=0; end
    if ind+1<=nall, a0(ind+1,ind)=0; end
end

return
%