function [it] = diffconv2d(eps,L)
% L is the number of discretization points in each direction (1,2)
I = speye(L,L);
E = sparse(2:L,1:L-1,1,L,L);
T = ???;
A = (kron(T,I)+kron(I,T))*eps*(L+1)^2 + ???;
b=ones([L*L,1]);
[sol,fl,rr,it]=gmres(A,b,1000,1.e-6,1000);
%plot sol
x=linspace(1/L,1-1/L,L);
y=linspace(1/L,1-1/L,L);
[xx, yy] = meshgrid (x, y);
k=0;
for i=1:L
 for j=1:L
  k=k+1;
  u(i,j)=sol(k);
  u(i,j)=sol(k);
 end
end
mesh(x,y,u);
end