function  uinfty=diffconv2dparab(eps,L,M)
% L is the number of discretization points in each direction (1,2)
% M is the number of time steps
tau=10/M;
I = speye(L,L);
E = sparse(2:L,1:L-1,1,L,L);
T = 2*I-E-E';
A = (kron(T,I)+kron(I,T))*eps*(L+1)^2+kron(I,E'-E)*(L+1)/2;
mat=speye(L*L,L*L)+???;
b=ones([L*L,1]);
u=zeros([L*L,1]);
for i=1:M;
  rhs=u+???;
  [lower,upper] = ilu (mat);
  u=gmres(mat,rhs,1000,1.e-6,1000,lower,upper,???);
end
%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;
  sol(i,j)=u(k);
  sol(i,j)=u(k);
 end
end
uinfty=max(abs(u));
mesh(x,y,sol);
end