function [u,uder] = intbvp_system(P,f,Wa,Wb,r)
% solves a system of differential equations in integral form
% solves a system of boundary value problems:
%
% P_{m}(y) u^{(m)}(y) + P_{m-1}(y) u^{(m-1)}(y) + ... + 
%   P_{0}(y) u(y) = f(y),
%
% where
%
% P is a matrix of functions
%
% a < y < b
%
% sum_{j=0}^{m-1} [Wa_{ij} u^{(j)}(a) + Wb_{ij} u^{(j)}(b)] = r_{i}
% i = 1, ... , mu
%
%
% written by: Binh Lieu
%
% reference: T.A. Driscoll / Journal of Computational Physics 229 (2010)
% pp. 5980 - 5998
%
% additional m-file needed: intdata_system.m
n = size(P,1);  % number of states in your system of ODE

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

dd = domain(P{1,1});     % domain of your functions

fzero = chebfun(0,dd);

[AT,ACm,WT,WCm,C,T,H] = intdata_system(P,Wa,Wb);

% Constructing operator form of ACm and H from cell format
Hn = [];
ACmn = [];

for ii = 1:n
    ACmjj = [];
    Hjj = [];
    for jj = 1:n
        ACmjj = [ACmjj, ACm{ii,jj}];
        Hjj = [Hjj, H{ii,jj}];
    end
    
    ACmn = [ACmn; ACmjj];
    Hn = [Hn; Hjj];
end

ACm = ACmn;
H = Hn;

% Construct chebfun solution for v = D^m * u, and integration constants.
for ii = 1:n
    bcr(:,ii) = fzero;
end
for jj = 1:n
    rr = WT{jj}\r{jj};
    for ii = 1:n
        ATij = AT{ii,jj};
        bcr(:,ii) = bcr(:,ii) + ATij*rr;
    end
end

vv = (ACm - H)\(f - bcr);

for ii = 1:n
    kk{ii} = WT{ii}\(r{ii} - WCm{ii}*vv(:,ii));
end

% Reconstruct all derivatives of the solution accurately.
for ii = 1:n
    Cju = vv(:,ii);
    ud = vv(:,ii);
    kki = kk{ii};
    Ti = T{ii};
    
    for jj = 1:m(ii)
        Cju = C*Cju;
        ud(:,jj + 1) = Cju + Ti(:,1:jj)*kki((m(ii)-jj+1):m(ii));
    end
    
    uder{ii} = ud;
    u{ii} = ud(:,end);
end

end

function [AT,ACm,WT,WCm,C,Tn,H] = intdata_system(P,Wa,Wb)
% ------------------------------------------------------------------------
% creates system's operators to be used in intbvp or inteigs
%
% P_{m}(y) u^{(m)}(y) + P_{m-1}(y) u^{(m-1)}(y) + ... + 
%   P_{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];
% P         -   a cell-array containing quasimatrices of functions;
% Wa, Wb    -   cell-arrays containing the boundary condition matrices;
%
% written by: Binh Lieu, 2011
% references: 
%   1)  T.A. Driscoll / Journal of Computational Physics 229 (2010) 
%       pp. 5980 - 5998
% ------------------------------------------------------------------------
%% System parameters
n = size(P,1);  % number of states in your system of ODE

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

dd = domain(P{1,1});    % domain of your functions
b = dd(end);            % end points of the boundaries
C = cumsum(dd);         % indefinite integration operator

%% Formulate quasimatrix for null space of D^m.
T = chebfun(1,dd);

for k = 2:max(m)
    T(:,k) = C*T(:,k-1);
end

for jj = 1:n
    Tn{jj} = T(:,1:m(jj));
end

%% Interior operator A applied to T, by Horner's method.
for ii = 1:n
    for jj = 1:n
        ATij = 0; DT = Tn{jj};
        Pij = P{ii,jj};
        for ind = 1:m(jj)
            ATij = ATij + diag(Pij(:,ind))*DT;
            DT = [ chebfun(0,dd), DT(:,1:m(jj)-1) ];
        end
        AT{ii,jj} = ATij;
    end
end

%% Boundary operator W applied to T.
for jj = 1:n
    TT = Tn{jj};
    WTj = Wa{jj} + Wb{jj}*toeplitz([1 zeros(1,m(jj)-1)], TT(b,:));
    WT{jj} = WTj;
end

%% Operators applied to C^m integration operator, again using Horner.
Eb = linop(@(n) [zeros(1,n-1) 1],@(u) feval(u,b), dd);
for ii = 1:n
    for jj = 1:n
        Pij = P{ii,jj};
        ACmij = diag(Pij(:,1));
        for ind = 1:m(jj)-1
            ACmij = ACmij*C + diag(Pij(:,ind+1));
        end
        ACmij = ACmij*C + diag(Pij(:,m(jj)+1));
        ACm{ii,jj} = ACmij;
    end
end

for jj = 1:n
    Wbj = Wb{jj};
    WCmj = Wbj(:,1)*Eb;
    for ind = 1:m(jj)-1
        WCmj = WCmj*C + Wbj(:,ind+1)*Eb;
    end
    WCmj = WCmj*C;
    WCm{jj} = WCmj;
end

%% Chebop form of ``outer products".
for ii = 1:n
    for jj = 1:n
        
        ATij = AT{ii,jj};
        WTj = WT{jj};
        WCmj = WCm{jj};
        
        mat = @(n) ATij(chebpts(n,dd),:)*(WTj\WCmj(n)); % n-by-n realization
        op = @(v) ATij*(WTj\(WCmj*v)); % functional expression
        
        H{ii,jj} = linop(mat,op,dd);
        
    end
end

end