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))];
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+R1'*(A1\(R1*(f-A*udd)));      % first subdomain solve
  U(2:m+1,2:m+1)=reshape(udd,m,m);
  mesh(x,y,U)
  xlabel('x');ylabel('y');
  pause
  udd=udd+R2'*(A2\(R2*(f-A*udd)));      % second subdomain solve
end
