m=15;                                   % number of gridpoints  
A=Laplacian(m,2);                       % five point Laplacian
f=zeros(m*m,1); f(m:m:end)=1;           % put bc into the rhs
u=A\f;                                  % solve by sparse Gaussian elimination
ai=9; bi=7;                             % alpha=(ai+1)*h, beta=(bi-1)*h
R1=[speye(m*ai) sparse(m*ai,m*(m-ai))]; % construct restriction matrices
R2=[sparse(m*(m-bi+1),m*(bi-1)) speye(m*(m-bi+1))];
mi=floor((ai+bi)/2);
R1t=R1;R1t(m*mi+1:end,:)=0;             % up to the middle included
R2t=R2;R2t(1:(mi-bi+1)*m,:)=0;         
A1=R1*A*R1';                            % construct subdomain matrices 
A2=R2*A*R2';
h=1/(m+1); x=0:h:1; y=x;                % mesh point vectors
udd=zeros(size(f));                     % initial guess
U=zeros(m+2); U(end,1:m+2)=1;           % for plotting
for n=0:10
  err(n+1)=max(max(abs(u-udd)));        % compute error 
  U(2:m+1,2:m+1)=reshape(udd,m,m);
  mesh(x,y,U)
  xlabel('x');ylabel('y');
  pause
  udd=udd+R1t'*(A1\(R1*(f-A*udd)))+R2t'*(A2\(R2*(f-A*udd)));
end
