function [u,uder] = intbvp(P,f,Wa,Wb,r)
% differential equations in integral form
% solves a m-th order linear boundary value problem:
%
% p_{m}(y) u^{(m)}(y) + p_{m-1}(y) u^{(m-1)}(y) + ... + 
%   p_{0}(y) u(y) = f(y),
%
% a < y < b
%
% sum_{j=0}^{m-1} [wa_{ij} u^{(j)}(a) + wb_{ij} u^{(j)}(b)] = r_{i}
% i = 1, ... , m
%
% input:
%
% P(1) = p_{0}
% P(2) = p_{1}
% ...
% P(m+1) = p_{m}
%
% Wa(i,j) = wa_{ij}
% Wb(i,j) = wb_{ij}
%
% written by: Binh Lieu
%
% reference: T.A. Driscoll / Journal of Computational Physics 229 (2010)
% pp. 5980 - 5998
%
% additional m-file needed: intdata.m

m = size(P,2)-1; % order of ODE

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

% Construct chebfun solution for v = D^m * u, and integration constants.
v = (ACm - H)\(f - AT*(WT\r));
% v = gmres(ACm -H, f - AT*(WT\r));
k = WT\(r - WCm*v);

% Reconstruct all derivatives of the solution accurately.
Cjv = v; uder = v;

for j = 1:m
    Cjv = C*Cjv;
    uder(:,j + 1) = Cjv + T(:,1:j)*k((m-j+1):m);
end

u = uder(:,end); % same as u = C^m * v + T * k

end


%% generate operators for integral boundary value problem
function [AT,ACm,WT,WCm,C,T,H] = intdata(P,Wa,Wb)
% differential equations in integral form
% solves a m-th order linear boundary value problem:
%
% p_{m}(y) u^{(m)}(y) + p_{m-1}(y) u^{(m-1)}(y) + ... + 
%   p_{0}(y) u(y) = f(y),
%
% a < y < b
%
% sum_{j=0}^{m-1} [wa_{ij} u^{(j)}(a) + wb_{ij} u^{(j)}(b)] = r_{i}
% i = 1, ... , m
%
% input:
%
% P(1) = p_{0}
% P(2) = p_{1}
% ...
% P(m+1) = p_{m}
%
% Wa(i,j) = wa_{ij}
% Wb(i,j) = wb_{ij}
%
% written by: Binh Lieu
% reference: T.A. Driscoll / Journal of Computational Physics 229 (2010)
% pp. 5980 - 5998
%
% addition m-file needed: intbvp.m

m = size(P,2) - 1; % order of ODE
d = domain(P(:,1));
a = d(1); b = d(end);
C = cumsum(d);

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

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

end

% Interior operator A applied to T, by Horner's method.
AT = 0; DT = T;

for j = 1:m
    AT = AT + diag(P(:,j))*DT;
    DT = [ chebfun(0,d), DT(:,1:m-1) ];
end

% Boundary operator W applied to T.
WT = Wa + Wb*toeplitz([1 zeros(1,m-1)], T(b,:));

% Operators applied to C^m integration operator, again using Horner.
ACm = diag(P(:,1));
Eb = linop(@(n) [zeros(1,n-1) 1],@(u) feval(u,b), d);
WCm = Wb(:,1)*Eb;

for j = 1:m-1
    ACm = ACm*C + diag(P(:,j + 1));
    WCm = WCm*C + Wb(:,j + 1)*Eb;
end

ACm = ACm*C + diag(P(:,m + 1));
WCm = WCm*C;

% Chebop form of ``outer products".
mat = @(n) AT(chebpts(n,d),:)*(WT\WCm(n)); % n-by-n realization
op = @(v) AT*(WT\(WCm*v)); % functional expression
H = linop(mat,op,d);

end