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=10; bi=6;                     % alpha=ai*h, beta=bi*h, 2h less than RAS!
R1=[speye(m*ai) sparse(m*ai,m*(m-ai))];
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;         
p1=4.4760; p2=p1;                % put optimized values 
% p1=7.0775; p2=3.8679;            % two sided optimized values  
A1=R1*A*R1';
A2=R2*A*R2';                     % change to Robin transmission conditions
h=1/(m+1);                       % mesh size scaling needed  
A1(end-m+1:end,end-m+1:end)=0.5*A1(end-m+1:end,end-m+1:end)+p1*h*speye(m);
A2(1:m,1:m)=0.5*A2(1:m,1:m)+p2*h*speye(m);
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,x,U)
  xlabel('x');ylabel('y');
  pause
  udd=udd+R1t'*(A1\(R1*(f-A*udd)))+R2t'*(A2\(R2*(f-A*udd)));
end
