function [A,Aeig] = brussel(m); 
% Brusselator model as described in Saad's paper. 

m = m/2; 
%
h = 1/(m+1); 
L = 0.51302; 
dx = 0.008; 
dy = 0.004;  
alpha = 2; 
beta = 5.45; 
tmp = (h*L)^2; 
c1 = dx/tmp; 
d1 = dy/tmp; 

%T = zeros(m); 
%for j=1:m,
%T(j,j) = -2; 
%end
%for j=1:m-1,
%    T(j,j+1) = 1;
%end
%for j=2:m,
%    T(j,j-1) = 1;
%end 
%AA(1:m,1:m) = c1*T + (beta-1)*eye(m); 
%AA(1:m,m+1:2*m) = alpha^2*eye(m);
%AA(m+1:2*m,1:m) = -beta*eye(m); 
%AA(m+1:2*m,m+1:2*m) = d1*T - alpha^2*eye(m); 
%eigofaa = eig(AA);
%hold on
%plot(real(eigofaa),imag(eigofaa),'y+')

nz = 0; 
for i=1:m, 
  nz = nz+1; 
  ia(nz) = i; 
  ja(nz) = i;
  a(nz) = -2*c1 + (beta-1);
end  

for i=1:m-1, 
  nz = nz+1; 
  ia(nz) = i; 
  ja(nz) = i+1;
  a(nz) = c1;
end 
for i=2:m, 
  nz = nz+1; 
  ia(nz) = i; 
  ja(nz) = i-1;
  a(nz) = c1;
end 

for i=m+1:2*m,
    nz = nz + 1; 
    ia(nz) = i; 
    ja(nz) = i; 
    a(nz) = -2*d1 - alpha^2; 
end 
for i=m+1:2*m-1,
    nz = nz + 1;
    ia(nz) = i;
    ja(nz) = i+1;
    a(nz) = d1;
end
for i=m+2:2*m,
    nz = nz + 1;
    ia(nz) = i;
    ja(nz) = i-1;
    a(nz) = d1;
end

for i=1:m,
    nz = nz + 1;
    ia(nz) = i;
    ja(nz) = i+m;
    a(nz) = alpha^2;
end

for i=m+1:2*m,
    nz = nz + 1;
    ia(nz) = i;
    ja(nz) = i-m;
    a(nz) = -beta;
end

A = sparse(ia,ja,a); 

%compute all eigenvalues 
for j=1:m,
   eigofT(j) = -2*(1- cos(pi*j*h) );
end
Aeig = zeros(2*m,1); 
for j=1:m,
    coeff(1) = 1;
    coeff(2) = alpha^2 - (beta - 1) - (d1+c1)*eigofT(j);
    coeff(3) = beta*alpha^2 + c1*d1*eigofT(j)^2 + d1*(beta-1)*eigofT(j) - ...
               alpha^2*c1*eigofT(j) - alpha^2*(beta-1);
    d = roots(coeff);
    Aeig(j) = d(1);
    Aeig(m+j) = d(2);
end;
