function [Sfun,Sval,A,B,C,Wa,Wb,A0s,B0s,C0s] = ...
    svdfr(A0,B0,C0,Wa0,Wb0,LRfuns,Nsigs)
% ========================================================================
% [Sfun,Sval] = svdfr(A0,B0,C0,Wa0,Wb0,LRfuns,Nsigs)
%  
% Given a two point boundary value representation of the frequency response
% operator
%  
%       { A0*phi = B0*d,
%   T:  { u = C0*phi,
%       { 0 = Wa0*phi(a) + Wb0*phi(b),
%  
% solve the eigenvalue problem
%  
%       T*Ts*Sfun = Sval*Sfun, or Ts*T*Sfun = Sval*Sfun,
%  
% where Ts is the adjoint of the frequency response operator T
%  
%        { A0s*psi = B0s*f,
%   Ts:  { g = C0s*psi,
%        { 0 = Wa0s*psi(a) + Wb0s*psi(b).
%  
% Inputs:
%  
%   LRfuns  = 1     --> solve for left singular functions: T*Ts
%                   --> determine spatial profile of the output
%   LRfuns  = 0     --> solve for right singular functions: Ts*T
%                   --> determine spatial profile of the input
%  
%   Nsigs --> number of desired singular values (default: Nsigs = 1)
%  
% Outputs:
%   
%   Sval    --> singular values of T arranged in descending order
%   Sfun    --> singular functions associated with Sval
%  
% written by: Binh Lieu
%
% reference:
%  
% B. K. Lieu and M. R. Jovanovic, "Computation of frequency responses
% for linear time-invariant PDEs on a compact interval", Journal of 
% Computational Physics, 2013, doi:10.1016/j.jcp.2013.05.010;
% also arXiv:1112.0579
%  
% ========================================================================

% assign default values
if nargin < 7
    Nsigs = 1;
end

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%% Chebfun's functions and operators
dom = domain(A0{1,1});    % domain of your functions
fone = chebfun(1,dom);   % one function
fzero = chebfun(0,dom);  % zero function

a = dom(1); % left-boundary
b = dom(2); % right-boundary

% identity and zero operators
II = eye(dom);
ZZ = zeros(dom);

% differential operator
D1 = diff(dom,1);

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%% System parameters
% number of fields in the RHS
n = size(A0,1);

% determine the highest order of each state appearing the equations
m = zeros(n,1);
for jj = 1:n
    m(jj) = size(A0{jj,jj},2) - 1;
end

mmax = max(m);  % highest derivative appearing in the LHS

% determine the number of input fields and the highest derivative
nd = size(B0,2);
md = size(B0{1,1},2) - 1;

% determine the number of output fields and the highest derivative
nu = size(C0,1);
mu = size(C0{1,1},2) - 1;

% storage for alpha_i (coefficients of A arranged in matrix form)
avec = fzero;  % storage vector of functions
for jj = 1:n*n
    avec(:,jj) = fzero;
end

alphai = cell(mmax+1, 1);
for jj = 1:mmax+1
    alphai{jj} = avec;
end

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%% Adjoint system
% ------------------------------------------------------------------------
% determine the adjoint operators A0s, B0s, and C0s
A0s = cell(n,n);
B0s = cell(nd,n);
C0s = cell(n,nu);
A21 = cell(n,n);

for indi = 1:n
    
    % ..........................................................
    % A0s
    for indj = 1:n
        indij = n*(indi - 1) + indj;
        
        % place holder for the ji component of A0s
        if m(indi) > m(indj)
            A0sji = 0*A0{indj,indi};
        else
            A0sji = 0*A0{indi,indj};
        end
        
        % ij-component of A0
        A0ij = A0{indi,indj};
        
        % Construct the operator A21 (zero coefficients)
        % will be used to create the block operator 
        % A = [A0, B0*B0s; A21, A0s] or
        % A = [A0s, C0s*C0; A21, A0]
        A21{indi,indj} = 0*A0ij;
        
        % if A0ij has all zero coefficients, continue to next component
        if norm(A0ij) == 0
            A0s{indj,indi} = A0sji;
            continue
        end
        
        % looping through each derivative
        for indni = 1:m(indj)+1
            
            pm1 = (-1)^(indni-1);
             
            for kk = 1:indni
                % binomial coefficient for product rule
                bc = binomial(indni-1,kk-1);
                % summing all terms coming in a_{ij,k}
                A0sji(:,kk) = A0sji(:,kk) + ...
                    pm1*bc*diff(conj(A0ij(:,indni)),indni - kk);
            end
            
            % ..........................................................
            % rearranging A in matrix form 
            % (will be used to determine adjoint boundary conditions)
            avec = alphai{indni};
            avec(:,indij) = A0ij(:,indni);
            
            alphai{indni} = avec;
        end
        
        A0s{indj,indi} = A0sji;
    end
    
    % ..........................................................
    % B0s
    for indj = 1:nd
        % place holder for the ji component of B0s
        B0sji = 0*B0{indi,indj};
        
        % ij-component of A0
        B0ij = B0{indi,indj};
        
        % if B0ij has all zero coefficients, continue to next component
        if norm(B0ij) == 0
            B0s{indj,indi} = B0sji;
            continue
        end
        
        % looping through each derivative
        for indni = 1:md+1
            
            pm1 = (-1)^(indni-1);
             
            for kk = 1:indni
                % binomial coefficient for product rule
                bc = binomial(indni-1,kk-1); 
                B0sji(:,kk) = B0sji(:,kk) + ...
                    pm1*bc*diff(conj(B0ij(:,indni)),indni - kk);
            end
            
        end
        
        B0s{indj,indi} = B0sji;
    end
    
    % ..........................................................
    % C0s
    for indj = 1:nu
        % place holder for the ij component of C0s
        C0sij = 0*C0{indj,indi};
        
        % ij-component of A0
        C0ji = C0{indj,indi};
        
        % if C0ji has all zero coefficients, continue to next component
        if norm(C0ji) == 0
            C0s{indi,indj} = C0sij;
            continue
        end
        
        % looping through each derivative
        for indni = 1:mu+1
            
            pm1 = (-1)^(indni-1);
             
            for kk = 1:indni
                % binomial coefficient for product rule
                bc = binomial(indni-1,kk-1); 
                C0sij(:,kk) = C0sij(:,kk) + ...
                    pm1*bc*diff(conj(C0ji(:,indni)),indni - kk);
            end
            
        end
        
        C0s{indi,indj} = C0sij;
    end
    
end

% ------------------------------------------------------------------------
% determine the adjoint boundary conditions
% formula is taken from: 
% M. Renardy & R. C. Rogers
% An Introduction to Partial Differential Equations
% Springer, 2004, pp. 163 - 165 
% 1) Find the A_{\alpha} matrix at the boundary points a and b
Ab = zeros(n*mmax,n*mmax);
Aa = zeros(n*mmax,n*mmax);

for indi = 1:mmax
    for indj = 1:mmax-(indi-1)
  
        Abij = zeros(n,n);
        Aaij = zeros(n,n);
        
        for kk = 1:(mmax-indj-indi+2)
            bc = binomial(kk+indi-2,indi-1);
            avec = diff(conj(alphai{kk+indj+indi-1}),kk-1);
            
            amatb = vec2mat(avec(b,:),n).';
            amata = vec2mat(avec(a,:),n).';
        
            Abij = Abij + (-1)^(kk+indi-2)*bc*amatb;
            Aaij = Aaij + (-1)^(kk+indi-2)*bc*amata;
        end
        
        Ab((indi-1)*n+1:indi*n, (indj-1)*n+1:indj*n) = Abij;
        Aa((indi-1)*n+1:indi*n, (indj-1)*n+1:indj*n) = Aaij;
    end
end

% 2) Find the basis for the space X
% 3) Apply A_alpha to the X
% 4) Form boundary conditions matrices
Wb = cell(2*n,1);
Wa = cell(2*n,1);

% loop over all variables
for indn = 1:n
    Wbn = Wb0{indn};
    Wan = Wa0{indn};
    
    nb = size(Wbn,1);
    na = size(Wan,1);
    
    [Ub,Sb,Vb] = svd(Wbn);
    [Ua,Sa,Va] = svd(Wan);
    
    % place holder for adjoint boundary conditions
    Wbsn = zeros(m(indn) - nb, m(indn));
    Wasn = zeros(m(indn) - nb, m(indn));
    
    % determine boundary conditions at right end for adjoint system
    for kk = 1:m(indn) - nb
        Xb = zeros(n,mmax);
        
        Xb(indn,1:m(indn)) = Vb(:,nb+kk).';
        Xb = vec2mat(Xb.',1);
        
        Yb = Ab*Xb;
        Yb = (Yb(indn:n:end).')./max(abs(Yb));
        
        Wbsn(kk,:) = Yb(1:m(indn));
    end
    
    % determine boundary conditions at left end for adjoint system
    for kk = 1:m(indn) - na
        Xa = zeros(n,mmax);
        
        Xa(indn,1:m(indn)) = Va(:,na+kk).';
        Xa = vec2mat(Xa.',1);
        
        Ya = Aa*Xa;
        Ya = (Ya(indn:n:end).')./max(abs(Ya));
        
        Wasn(kk,:) = Ya(1:m(indn));
    end
    
    % Pad boundary conditions with zero to put them in 
    % required standard form
    % Boundary condition matrices for original and adjoint systems
    if LRfuns
        Wb{indn} = [zeros(na, m(indn)); Wbn];
        Wa{indn} = [Wan; zeros(nb, m(indn))];
        
        Wb{n+indn} = [zeros(nb, m(indn)); Wbsn];
        Wa{n+indn} = [Wasn; zeros(na, m(indn))];
    elseif ~LRfuns 
        Wb{indn} = [zeros(nb, m(indn)); Wbsn];
        Wa{indn} = [Wasn; zeros(na, m(indn))];
        
        Wb{n+indn} = [zeros(na, m(indn)); Wbn];
        Wa{n+indn} = [Wan; zeros(nb, m(indn))];
    else
        error('Invalid value of LRfuns. Please choose: \n LRfuns = 1 for left singular functions \n LRfuns = 0 for for right singular functions \n');
    end
end

% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
%% System coefficients for the composition frequency response operator 
% LRfuns = 1 --> left singular functions: T*Ts
% LRfuns = 0 --> right singular functions: Ts*T

% ------------------------------------------------------------------------
% find the operator T*Ts and determine the left singular functions
if LRfuns
    BBs = cell(n,n);
    % ........................................
    % determine the operator B0*B0s
    % loop over all the elements of B0*B0s, which has size n x n
    for indi = 1:n
        for indj = 1:n
            % place holder for ij component of B0*B0s
            BBsij = 0*A0s{1,indj};
            
            % loop through each component of the forcing
            % need to do a summation (matrix multiplication of B0*B0s)
            % need to take care of product rule when composing the
            % two operators
            for indni = 1:nd
                % ij component of B0
                B0ij = B0{indi,indni};
                
                % ij component of B0s
                B0sij = B0s{indni,indj};
                
                if norm(B0sij) == 0
                    continue
                end
                
                for indmi = 1:md+1
                    
                    for kk = indmi:-1:1
                        % binomial coefficient for product rule
                        bc = binomial(indmi-1,kk-1);
                        
                        BBsij(:,indmi-kk+1:indmi-kk+md+1) = ...
                            BBsij(:,indmi-kk+1:indmi-kk+md+1) + ...
                            bc*diag(B0ij(:,indmi))*diff(B0sij(:,:),kk-1);
                        
                    end
                    
                end
                
            end
            
            % we need to bring the BBs operator to LHS, hence the minus
            BBs{indi,indj} = -BBsij;
        end
    end
    
    % ........................................
    % determine the operator A, B, and C
    A = [A0, BBs; A21, A0s];
    
    C = cell(nu,1);
    B = cell(1,nu);
    
    % create a dummy variable with zero operators for padding
    Zrvec = ZZ;
    Zcvec = ZZ;
    
    for indj = 1:n-1
        Zrvec = [Zrvec, ZZ];
        Zcvec = [Zcvec; ZZ];
    end
    
    % determine operators B and C from C0s and C0
    for indi = 1:nu
        Ci = [];
        Bi = [];
        
        for indj = 1:n
            C0ij = C0{indi,indj};
            C0sji = C0s{indj,indi};
            
            % initiate component of Cij with first coefficient of C0ij
            Cni = diag(C0ij(:,1));
            Bni = diag(C0sji(:,1));
            
            for indni = 1:mu
                Cni = Cni + diag(C0ij(:,indni+1))*(D1^(indni));
                Bni = Bni + diag(C0sji(:,indni+1))*(D1^(indni));
            end
            
            Ci = [Ci, Cni];
            Bi = [Bi; Bni];
        
        end
        
        % pad end components with zeros
        C{indi,1} = [Ci, Zrvec];
        B{1,indi} = [Zcvec; Bi];
    end
    
% ------------------------------------------------------------------------
% find the operator Ts*T and determine the right singular functions
elseif ~LRfuns
    CsC = cell(n,n);
    % ........................................
    % determine the operator C0^{\star}*C0
    % loop over all the elements of C0^{\star}*C0, which has size n x n
    for indi = 1:n
        for indj = 1:n
            % place holder for ij component of C0^{\star}*C0
            CsCij = 0*A0s{1,indj};
            
            % loop through each component of the forcing
            % need to do a summation (matrix multiplication of
            % C0^{\star}*C0)
            % need to take care of product rule when composing the
            % two operators
            for indni = 1:nu
                % ij component of C0^{\star}
                C0sij = C0s{indi,indni};
                
                % ij component of C0s
                C0ij = C0{indni,indj};
                
                if norm(C0ij) == 0
                    continue
                end
                
                for indmi = 1:mu+1
                    
                    for kk = indmi:-1:1
                        % binomial coefficient for product rule
                        bc = binomial(indmi-1,kk-1);
                        
                        CsCij(:,indmi-kk+1:indmi-kk+md+1) = ...
                            CsCij(:,indmi-kk+1:indmi-kk+md+1) + ...
                            bc*diag(C0sij(:,indmi))*diff(C0ij(:,:),kk-1);
                        
                    end
                    
                end
                
            end
            
            % need to bring CsC operator to LHS, hence the minus
            CsC{indi,indj} = -CsCij;
        end
    end
    
    % ........................................
    % determine the operator A, B, and C
    A = [A0s, CsC; A21, A0];
    
    C = cell(nd,1);
    B = cell(1,nd);
    
    % create a dummy variable with zero operators for padding
    Zrvec = ZZ;
    Zcvec = ZZ;
    
    for indj = 1:n
        Zrvec = [Zrvec, ZZ];
        Zcvec = [Zcvec; ZZ];
    end
    
    % determine operators B and C from B0 and B0s
    for indi = 1:nd
        Ci = [];
        Bi = [];
        
        for indj = 1:n
            B0ij = B0{indi,indj};
            B0sji = B0s{indj,indi};
            
            % initiate component of Cij with first coefficient of B0sji
            % initiate component of Bij with first coefficient of B0ij
            Cni = diag(B0sji(:,1));
            Bni = diag(B0ij(:,1));
            
            for indni = 1:md
                Cni = Cni + diag(B0sji)*(D1^(indni));
                Bni = Bni + diag(B0ij)*(D1^(indni));
            end
            
            Ci = [Ci, Cni];
            Bi = [Bi; Bni];
        
        end
        
        % pad the end components with zeros
        C{indi,1} = [Ci, Zrvec];
        B{1,indi} = [Zcvec; Bi];
    end
    
else
    error('Invalid value for LRfuns. Please choose: \n LRfuns = 1 for left singular functions \n LRfuns = 0 for for right singular functions \n');
end

%% Computing singular values
[xint,lam] = inteigs_system_TTs(A,Wa,Wb,B,C,Nsigs,'LR');

Sval = sqrt(real(lam));
Sfun = xint;

end

function [v,lam] = inteigs_system_TTs(A,Wa,Wb,B,C,Neigs,sigma)
% ========================================================================
% [v,lam] = inteigs_system_TTs(A,Wa,Wb,B,C,Neigs,sigma)
%  
% Given a system of differential equations,
% Solve the eigenvalue problem in integral form:
%  
%       T*v = lam*v,
%  
% where T is an operator given by
%  
%       { A*u = B*f,
%   T:  { v = C*u,
%       { 0 = Wa*u(a) + Wb*u(b),
%  
% where
%  
% A_{m}(y) u^{(m)}(y) + A_{m-1}(y) u^{(m-1)}(y) + ... + 
%   A_{0}(y) u(y) = B*f, a < y < b
%  
% sum_{j=0}^{m-1} [Wa_{ij} u^{(j)}(a) + Wb_{ij} u^{(j)}(b)] = 0, 
%   i = 1, ... , mu
%  
% Inputs:
%  
%   Neigs   --> number of desired eigenvalues (default: Neigs = 6)
%   sigma   --> specifies criterion for determining eigenvalues 
%               (default: sigma = 'LR')
%               'LR' - eigenvalues with largest real parts
%  
% Outputs:
%   
%   lam     --> eigenvalues of T
%   v       --> eigenfunctions associated with lam
%  
% written by: Binh Lieu, 2011
% references:
%  
%   1)  B.K. Lieu & M.R. Jovanovic, "Computation of frequency responses of 
%       linear time-invariant PDEs on a compact interval", 
%       Journal of Computational Physics, submitted (2011).
%  
%   2)  T.A. Driscoll, "Automatic spectral collocation for integral, 
%       integro-differential, and integrally reformulated differential 
%       equations", Journal of Computational Physics, 
%       229, pp. 5980 - 5998, (2010).
%  
% ========================================================================

if nargin < 6
    Neigs = 6;
    sigma = 'LR';
elseif nargin < 7
    sigma = 'LR';
end

%% System parameters
n = size(A,1);

% determine the highest order of each state appearing the equations
m = zeros(n,1);
for jj = 1:n
    m(jj) = size(A{jj,jj},2) - 1;
end

dd = domain(A{1,1});    % domain of your functions
fone = chebfun(1,dd);   % one function

%% Identity functional operator: Qu = u
% Need this to compute RHS of eigenvalue problem
Q = cell(n,n);
for ii = 1:n
    for jj = 1:n
        Qij = 0*A{ii,jj};
        
        if ii == jj
            Qij(:,1) = fone;
        end
        
        Q{ii,jj} = Qij;
    end
end
    
%% LHS of eigenvalue problem
[AK,AJm,WK,WJm,J,Kn,H,WKsingular] = intdata_system(A,Wa,Wb);

% if WK is not singular, continue with method 1
% this method can find an explicit expression for the integration constants
% in terms of the new variable from the boundary conditions
% if WK is singular, then we need to use an alternative method
% all necessary operators are computed in intdata_system
if ~WKsingular
    % Constructing operators form of AJm and H from cell format
    Hn = [];
    AJmn = [];
    
    for ii = 1:n
        AJmjj = [];
        Hjj = [];
        
        for jj = 1:n
            AJmjj = [AJmjj, AJm{ii,jj}];
            Hjj = [Hjj, H{ii,jj}];
        end
        
        AJmn = [AJmn; AJmjj];
        Hn = [Hn; Hjj];
    end
    
    AJm = AJmn;
    H = Hn;
    
    %% Determine operators on the RHS of eigenvalue problem
    % Note that in this problem, RHS is identity
    % However, we need the operators in the correct format
    [QK,QJm,WK,WJm,J,Kn,G] = intdata_system(Q,Wa,Wb);
    
    % Constructing operator form of QCm and G from cell format
    Gn = [];
    QJmn = [];
    
    for ii = 1:n
        QJmjj = [];
        Gjj = [];
        
        for jj = 1:n
            QJmjj = [QJmjj, QJm{ii,jj}];
            Gjj = [Gjj, G{ii,jj}];
        end
        
        QJmn = [QJmn; QJmjj];
        Gn = [Gn; Gjj];
    end
    
    QJm = QJmn;
    G = Gn;
    
    %% Construct chebfun operators for TTs
    Abar = AJm - H;
    Cbar = QJm - G;
    
    ninput = size(B,2);
    noutput = size(C,1);
    
    TTs = [];
    
    for ii = 1:noutput
        
        TTsj = [];
        Ci = C{ii,1};
        
        for jj = 1:ninput
            
            Bj = B{1,jj};
            
            % n-by-n realization
            mat = @(nn) Ci(nn)*Cbar(nn)*(Abar(nn)\Bj(nn));
            % functional expression
            op = @(vv) Ci*Cbar*(Abar\Bj)*vv;
            
            TTsj = [TTsj, linop(mat,op,dd)];
            
        end
        
        TTs = [TTs; TTsj];
        
    end
    
else
    % WK is singular, use method 2 where H is the operator containing all
    % information
    ninput = size(B,2);
    noutput = size(C,1);
    
    TTs = [];
    
    for ii = 1:noutput
        
        TTsj = [];
        Ci = C{ii,1};
        
        for jj = 1:ninput
            
            Bj = B{1,jj};
            
            % n-by-n realization
            mat = @(nn) Ci(nn)*H(nn)*Bj(nn);
            % functional expression
            op = @(vv) Ci*H*Bj*vv;
            
            TTsj = [TTsj, linop(mat,op,dd)];
            
        end
        
        TTs = [TTs; TTsj];
        
    end
    
end



%% Finding eigenvalues and eigenfunctions of T
II = eye(dd);
Ibig = blkdiag(II,ninput);

[v,L] = eigs(TTs,Ibig,Neigs,sigma);

lam = diag(L);

end

%% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% Subfunction to help set up system's operators
function [AK,AJm,WK,WJm,J,Kn,H,WKsingular] = intdata_system(A,Wa,Wb)
% ========================================================================
% Creates system's operators to be used in inteigs_system_TTs
%
% A_{m}(y) u^{(m)}(y) + A_{m-1}(y) u^{(m-1)}(y) + ... + 
%   A_{0}(y) u(y) = f(y),
%
% sum_{j=0}^{m-1} [Wa_{ij} u^{(j)}(a) + Wb_{ij} u^{(j)}(b)] = 0
% i = {1, ... , mu},
%
% where
%
% y         -   spatial variable with a < y < b;
% u         -   vector of state variables: u = [u_1; u_2, ... ; u_n];
% A         -   a cell-array containing quasimatrices of functions;
% Wa, Wb    -   cell-arrays containing the boundary condition matrices;
%
% written by: Binh Lieu, 2011
% references:
%  
%   [1] B.K. Lieu & M.R. Jovanovic, "Computation of frequency responses of 
%       linear time-invariant PDEs on a compact interval", 
%       Journal of Computational Physics, submitted (2011).
%  
%   [2] T.A. Driscoll, "Automatic spectral collocation for integral, 
%       integro-differential, and integrally reformulated differential 
%       equations", Journal of Computational Physics, 
%       229, pp. 5980 - 5998, (2010).
% ========================================================================
%% System parameters
n = size(A,1);  % number of states in your system of ODEs

% determine the highest order of each state appearing in the equations
m = zeros(n,1);
for jj = 1:n
    m(jj) = size(A{jj,jj},2) - 1;
end

dom = domain(A{1,1});    % domain of your functions
b = dom(end);            % end points of the left boundary
J = cumsum(dom);         % indefinite integration operator

I = eye(dom);
In = blkdiag(I,n);

%% Formulate quasimatrix for null space of D^m.
K = chebfun(1,dom);
Kn = cell(n,n);

for k = 2:max(m)
    K(:,k) = J*K(:,k-1);
end

for jj = 1:n
    for ii = 1:n
        if ii == jj
            Kn{ii,jj} = K(:,1:m(jj));
        else
            Kn{ii,jj} = 0*K(:,1:m(jj));
        end
    end
end

%% Interior operator A applies to K.
% Same operator as L12 in our paper [1]
AK = cell(n,n);
for ii = 1:n
    for jj = 1:n
        AKij = 0; DK = Kn{jj,jj};
        Aij = A{ii,jj};
        for ind = 1:m(jj)
            AKij = AKij + diag(Aij(:,ind))*DK;
            DK = [ chebfun(0,dom), DK(:,1:m(jj)-1) ];
        end
        AK{ii,jj} = AKij;
    end
end

%% Boundary operator W applies to K.
% same operator as L22 in our paper [1]
WK = cell(n,1);
WKsingular = 0;
for jj = 1:n
    KK = Kn{jj,jj};
    WKj = Wa{jj} + Wb{jj}*toeplitz([1 zeros(1,m(jj)-1)], KK(b,:));
    WK{jj} = WKj;
    
    if rank(WKj) < m(jj)
        WKsingular = 1;
    end
end

%% Operators apply to J^m integration operator.
% Operator A applies to the integration operator J, same operator as L11 in
% our paper [1]
Eb = linop(@(n) [zeros(1,n-1) 1],@(u) feval(u,b), dom);
AJm = cell(n,n);
for ii = 1:n
    for jj = 1:n
        Aij = A{ii,jj};
        AJmij = diag(Aij(:,1));
        for ind = 1:m(jj)-1
            AJmij = AJmij*J + diag(Aij(:,ind+1));
        end
        AJmij = AJmij*J + diag(Aij(:,m(jj)+1));
        AJm{ii,jj} = AJmij;
    end
end

% Boundary operator W applies to the integration operator J, 
% same operator as L21 in our paper [1]
WJm = cell(n,n);
for jj = 1:n
    Wbj = Wb{jj};
    WJmj = Wbj(:,1)*Eb;
    for ind = 1:m(jj)-1
        WJmj = WJmj*J + Wbj(:,ind+1)*Eb;
    end
    WJmj = WJmj*J;
    for ii = 1:n
        if ii == jj
            WJm{jj,ii} = WJmj;
        else
            WJm{jj,ii} = zeros(m(jj),m(jj))*WJmj;
        end
    end
end


%% Chebop form of ``outer products".
if WKsingular ~= 1
    % Matrix WK is not singular, hence, we can proceed with the original
    % method.
    H = cell(n,n);
    for ii = 1:n
        for jj = 1:n
            
            AKij = AK{ii,jj};
            WKj = WK{jj};
            WJmj = WJm{jj,jj};
            
            mat = @(N) AKij(chebpts(N,dom),:)*(WKj\WJmj(N)); % n-by-n realization
            op = @(v) AKij*(WKj\(WJmj*v)); % functional expression
            
            H{ii,jj} = linop(mat,op,dom);
            
        end
    end
    
else
    % Matrix WK is singular, we cannot proceed with the original method.
    % We cannot determine the constant of integration independently. Will
    % need to use another method that requires a different operator.
    
    % Build the operator AJ from cell format
    AJmn = [];
    for ii = 1:n
        AJmjj = [];
        
        for jj = 1:n
            AJmjj = [AJmjj, AJm{ii,jj}];
        end
        
        AJmn = [AJmn; AJmjj];
    end
    AJm = AJmn;
    
    % Build the operator AJ\AK and WJ*AJ\AK 
    Gm = cell(n,n);
    for jj = 1:n
        
        AKj = [AK{:,jj}];
        
        % Get the jj column of the operator AJ\AK
        mj = m(jj);
        AJAKj = 0*AKj;
        for indm = 1:mj;
            AJAKj(:,(indm-1)*n+1:indm*n) = (AJm\AKj(:,indm:mj:end));
        end
        
        for ii = 1:n
            WJmi = WJm{ii,ii};
            
            WJAJAKij = WJmi*AJAKj(:,ii:n:end);
            
            % Store WJAJAK temporarily in Gm
            Gm{ii,jj} = -WJAJAKij;
        end
    end
    
    % Get the matrix G = WK - WJ*AJ\AK
    WJmj = [WJm{1,:}];
    Gm{1,1} = WK{1} + Gm{1,1};
    G = [Gm{1,:}]; % Make Gm a matrix from cell array
    
    % Form the operator WJ*AJ\I
    mat = @(n) WJmj(n)*(AJm(n)\In(n)); % n-by-n realization
    op = @(v) WJmj*(AJm\(In*v)); % functional expression
    WJAJI = linop(mat,op,dom);
        
    for jj = 2:n
        Gm{jj,jj} = WK{jj} + Gm{jj,jj};
        G = [G; Gm{jj,:}];
        
        WJmj = [WJm{jj,:}];
        
        % Form the operator WJ*AJ\I
        mat = @(n) WJmj(n)*(AJm(n)\In(n)); % n-by-n realization
        op = @(v) WJmj*(AJm\(In*v)); % functional expression
        WJAJI = [WJAJI; linop(mat,op,dom)];
    end
    
    if rank(G) < sum(m)
        error('System is not well-posed. Please double check your system.');
    end

    % Build the operator AK*(G\WJ)*(AJ\I) and K*(G\WJ)*(AJ\I)
    AKi = [AK{1,:}];

    mat = @(n) AKi(chebpts(n),:)*(G\WJAJI(n)); % n-by-n realization
    op = @(v) AKi*(G\WJAJI)*v; % functional expression
    AKGWJAJI = linop(mat,op,dom);
    
    Kni = [Kn{1,:}];

    mat = @(n) Kni(chebpts(n),:)*(G\WJAJI(n)); % n-by-n realization
    op = @(v) Kni*(G\WJAJI)*v; % functional expression
    KGWJAJI = linop(mat,op,dom);
    
    Jbig = J^(m(1));
    
    for ii = 2:n
        AKi = [AK{ii,:}];
        Kni = [Kn{ii,:}];
        
        % Form the operator AK*(G\WJ)*(AJ\I)
        mat = @(n) AKi(chebpts(n),:)*(G\WJAJI(n)); % n-by-n realization
        op = @(v) AKi*(G\WJAJI)*v; % functional expression
        AKGWJAJI = [AKGWJAJI; linop(mat,op,dom)];
        
        % Form the operator K*(G\WJ)*(AJ\I)
        mat = @(n) Kni(chebpts(n),:)*(G\WJAJI(n)); % n-by-n realization
        op = @(v) Kni*(G\WJAJI)*v; % functional expression
        KGWJAJI = [KGWJAJI; linop(mat,op,dom)];
        
        Jbig = blkdiag(Jbig,J^(m(ii)));
    end
    
    % Build the operator H1 = Jbig*(AJ\(I + AK*(G\WJ)*(AJ\I)))
    mat = @(n) Jbig(n)*(AJm(n)\(In(n) + AKGWJAJI(n))); % n-by-n realization
    op = @(v) Jbig*(AJm\(In + AKGWJAJI))*v; % functional expression
    H1 = linop(mat,op,dom);
    
    % Build the operator H2 = K*(G\WJ)*(AJ\I)
    mat = @(n) In(n)*KGWJAJI(n); % n-by-n realization
    op = @(v) In*KGWJAJI*v; % functional expression
    H2 = linop(mat,op,dom);
    
    % Build the operator H = H1 + H2
    mat = @(n) H1(n) - H2(n); % n-by-n realization
    op = @(v) (H1 - H2)*v; % functional expression
    H = linop(mat,op,dom);
    
end
    
    
end


%% ========================================================================
% binomial.m 
% Written by David Terr, Raytheon, 5-11-04
% Given nonnegative integers n and m with m<=n, compute the 
% binomial coefficient n choose m. 
function bc = binomial(n,m)

if 2*m <= n
    m1 = m;
else
    m1 = n-m;
end

if n <= 100
    size = abs( n-m1 ) + 1;
    ps = pascal( size );
    bc = ps( size, m1+1 );
else    % suggested by Greg von Winckel to improve speed for large input
    bc = exp( gammaln(n+1) - gammaln(n-m+1) - gammaln(m+1) );
end

end


%% ========================================================================
function [mat, padded] = vec2mat(vec, matCol, padding)
%VEC2MAT Convert a vector into a matrix.
%   MAT = VEC2MAT(VEC, MATCOL) converts the vector VEC into a matrix with
%   MATCOL columns, creating one row at a time. If the length of VEC is not
%   a multiple of MATCOL, then the function places extra entries of 0 in the
%   last row of MAT.
%
%   MAT = VEC2MAT(VEC, MATCOL, PADDING) is the same as the first syntax,
%   except that the extra entries are taken from the matrix PADDING, in order.
%   If the number of elements in PADDING is insufficient, the last element is
%   used for the remaining entries.
%
%   [MAT, PADDED] = VEC2MAT(...) returns an integer PADDED that indicates
%   how many extra entries were placed in the last row of MAT.
%
%   See also RESHAPE.

%   Copyright 1996-2008 The MathWorks, Inc.
%   $Revision: 1.10.4.3 $

error(nargchk(2, 3, nargin,'struct'));	% 2 or 3 inputs required
if ndims(vec) > 2
    error('comm:vec2mat:InvalidVec','VEC cannot be an ND array.');
elseif (length(matCol) ~= 1 || ~isfinite(matCol) || ~isreal(matCol)...
        || floor(matCol) ~= matCol || matCol < 1)
    error('comm:vec2mat:InvalidMatcol','MATCOL must be a positive integer.');
end

[vecRow, vecCol] = size(vec);
vecLen = vecRow*vecCol;
if vecCol == matCol
    mat = vec;
    padded = 0;
    return;			% nothing to do
elseif vecRow > 1
    vec = reshape(vec.', 1, vecLen);
end

try
    if nargin < 3 || isempty(padding)
        padding = cast(0, class(vec));	% default padding
    else
        padding = cast(padding(:).', class(vec));
    end
catch exception
    throw(exception)
end
paddingLen = length(padding);

matRow = ceil(vecLen/matCol);
padded = matRow*matCol - vecLen;	% number of elements to be padded
vec = [vec, padding(1:min(padded, paddingLen)),...
       repmat(padding(paddingLen), 1, padded-paddingLen)];	% padding
mat = reshape(vec, matCol, matRow).';

end