function [yfe,ypfe,yppfe] = solve_1D_analytical(xgrid,xknots,sourcefunction,sporder,ngauss,BCrho,BCval)

% Spline matrix
[SS,SSp,SSpp] = splineval(xgrid,xknots,sporder,[1 0 0]);

nsp = size(SS,1);
[xgauss,~] = GaussGrid(xknots,ngauss);

MM.fixed = ones(size(xgauss));
DD.fixed = ones(size(xgauss));
KK.fixed = ones(size(xgauss));

% Mass matrix
di=0; dj=0; % zeroth derivative: MASS matrix
[MM.XI,MM.XIindex,MM.nsp] = FEconstruct(xknots,MM.fixed,sporder,ngauss,di,dj,[1 0 0]);

di=1; dj=1; % first derivative: STIFFNESS matrix
[DD.XI,DD.XIindex,DD.nsp] = FEconstruct(xknots,DD.fixed,sporder,ngauss,di,dj,[1 0 0]);

di=0; dj=[]; % for SOURCE term;
[KK.XI,KK.XIindex,KK.nsp] = FEconstruct(xknots,KK.fixed,sporder,ngauss,di,dj,[1 0 0]);

%% Evaluate FE matrices using matrix multiplication
% variable part of integrand on xgauss points
mgauss = ones(size(xgauss));
dgauss = -xgauss;
sgauss= sourcefunction(xgauss);

boundaryfunction = @(x) -10*(x-1);
bgauss = -boundaryfunction(xgauss);

%% allocate
M = zeros(nsp,nsp);
D = zeros(nsp,nsp);
S = zeros(nsp,1);
B = zeros(nsp,1);
%% compute matrices
M(MM.XIindex) =  MM.XI*mgauss;
D(DD.XIindex) =  DD.XI*dgauss;
S(KK.XIindex) =  KK.XI*sgauss;

B(KK.XIindex) =  KK.XI*bgauss;

%% Add boundary condition


nbc = numel(BCval);

BClhs = splineval(BCrho,xknots,sporder,[1 0 0])';
DBC = [D(1:end-nbc,1:end);-BClhs];
rhs = [S(1:end-nbc);BCval];

coef = -DBC\rhs;

yfe = SS'*coef ;
ypfe = SSp'*coef;
yppfe = SSpp'*coef;

end