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

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; 
A2G = [-1 ; zeros(n-1,1)]/h^2; 
AG1 = A1G'; 
AG2 = A2G';
AGG = 2/h^2; 

% Solve for the glue: 
A11f1 = KSolve(f1);         % Solve A11 \ f1
A22f2 = KSolve(f2);         % Solve A22 \ f2
A11A1G = KSolve(A1G);       % Solve A11 \ A1G 
A22A2G = KSolve(A2G);       % Solve A22 \ A2G
b = fG -  AG1 * (A11f1) - AG2 * (A22f2);    
uG = ( AGG - AG1 * (A11A1G) - AG2 * (A22A2G) ) \ b; % Solve for the 'glue'

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

% Final solution. Do not forget to pad with zeros (for the bcs) 
u = real([0 ; u1 ; uG ; u2 ; 0]);

end

function x = KSolve( 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