function u = domaindecomposition( n, rhs ) 
% Solve -u_xx = f on [0 2] with zero Dirichlet conditions by domain 
% decomposition. 
% 
% Alex Townsend, September 2014, updated September 2015.

h = 1/(n+1);    % grid spacing
x = linspace(0, 2, 2*n+3)'; x(1)=[]; x(end)=[];     % internal grid nodes.
f = rhs(x); f1 = f(1:n); fG = f(n+1); f2 = f(n+2:2*n+1); % righthand side.
  
% Get components that glue the system together. 
A1G = [zeros(n-1,1) ; -1]/h^2;   AG1 = A1G'; 
A2G = [-1 ; zeros(n-1,1)]/h^2;   AG2 = A2G';
AGG = 2/h^2; 

% Solve for the glue, schur complement: 
uG = ( AGG - AG1 * K(A1G) - AG2 * K(A2G) ) \ ... 
                             ( fG -  AG1 * K(f1) - AG2 * K(f2) );

% Solve the two problems independently, with DSTs:
u1 = K(f1 - A1G*uG); % Solve left domain
u2 = K(f2 - A2G*uG); % Solve right domain

% Final solution:
u = real([0 ; u1 ; uG ; u2 ; 0]); % pad with zeros for the bcs 

end

function x = K( b )
% Solve Kx = b, where K = 2nd order FD matrix. 
n = size(b, 1); h = 1/(n+1); 
s = sqrt(2/(n+1));                     % dst normalization
d = (2-2*cos((1:n)'*pi/(n+1)));        % Eigenvalues of K
x = s*dst1( s*dst1( b )./d )*h^2;      % Use eigenvalue decomposition
end 
