m=15;                                % number of gridpoints
h=1/(m+1);                           % include the h^2 scaling for
A=Laplacian(m,2)/h^2;                % the true five point Laplacian
f=zeros(m*m,1); f(m:m:end)=1/h^2;    % include h^2 scaling also in rhs
u=A\f;                               % since true solver Solve2dR is used
a=6;                                 % alpha=ai*h
F1=zeros(m,a+1); F1(m,:)=1/h^2;      % f for subdomain solves must
F2=zeros(m,m-a+2); F2(m,:)=1/h^2;    % include Robin bc columns
x1=0:h:a*h; x2=a*h:h:1; y=0:h:1;     % finite difference meshes
z1=zeros(1,a+1);z2=zeros(1,m-a+2);   % for plotting purposes
o1=ones(1,a+1);o2=ones(1,m-a+2);     % for plotting purposes
e=ones(m,1);                         % construct normal derivative
Na=[speye(m) -spdiags([-e +4*e -e]/2,[-1 0 1],m,m)]/h;
Nb=[-spdiags([-e +4*e -e]/2,[-1 0 1],m,m) speye(m)]/h;
la=zeros(m,1);                       % zero initial guess
U2=zeros(size(F2));
f1=@(al,L) (tanh(pi/L*al)+tanh(pi/L*(1-al)))*(coth(pi/L*al)+coth(pi/L*(1-al)));
th=2/(4+f1(a*h,1));                  % optimized relaxation parameter
pe=1e10;                             % 1e10 to emulate a Dirichlet condition
g=zeros(m,1);                        % left and right Dirichlet condition
for n=0:10
  err(n+1)=norm(la-u((a-1)*m+1:a*m),'inf');    
  U1=Solve2dR(F1,0,a,g*pe,la*pe,pe,pe);
  U2=Solve2dR(F2,a,m+1,la*pe,g*pe,pe,pe);
  mesh(x1,y,[z1;U1;o1]); hold on;mesh(x2,y,[z2;U2;o2]);hold off
  xlabel('x'); ylabel('y'); zlabel('Neumann-Neumann Iterate');
  axis([0 1 0 1 0 1])
  pause
  tb=Nb*[U2(:,1);U2(:,2)]+F2(:,1)*h/2;  
  ta=Na*[U1(:,end-1);U1(:,end)]+F1(:,end)*h/2;  
  psi1=Solve2dR(zeros(size(F1)),0,a,pe*g,ta+tb,pe,0);
  psi2=Solve2dR(zeros(size(F2)),a,m+1,ta+tb,g*pe,0,pe);
  la=la+th*(psi1(:,end)+psi2(:,1));
  mesh(x1,y,[z1;psi1;z1]); hold on;mesh(x2,y,[z2;psi2;z2]);hold off
  xlabel('x'); ylabel('y'); zlabel('Neumann-Neumann Correction');
  pause
end
