function error = IGAcoll_laplace1d(p,ncp,a,b,c)

% A (very simple) 1D IGA collocation code
%
% Solve a*u'' + b*u' + c*u = f with homogeneous Dirichlet b.c.'s and 
% with a manufactured body load s.t. the analytical solution is available
% 
% Input:  p    = approximation degree
%         ncp  = number of collocation points
%
% Output: error = (approximation of the) L2-norm error
%
%
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         % By:                                            %
%         %          Alessandro Reali                      %
%         %          University of Pavia                   %
%         % email:   alessandro.reali@unipv.it             %
%         % webpage: http://www.unipv.it/alereali          %
%         %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% set mesh [linear parameterization, continuity C^(p-1)]

knot      = augknt([0 1],p+1);
h         = 1/(ncp - p);
[~,knot] = bspkntins(p,0:1/p:1,knot,h:h:1-h);

% compute collocation points (Greville abscissae)
collpt = aveknt(knot,p+1);

% "loop" over active collocation points
ii = 2:ncp-1;
% set active collocation points
collpt = collpt(ii);
% compute basis functions and their derivatives at collocation pts
N = spcol(knot,p+1,sort([collpt,collpt,collpt]));
% construct the stiffness matrix
k_gl =  a*N(3:3:end,ii) + b*N(2:3:end,ii) + c*N(1:3:end,ii);
% construct the r.h.s.
f_gl = (c - a*4*pi^2)*sin(2*pi*collpt)' + b*2*pi*cos(2*pi*collpt)';

% solve linear system
sol = k_gl\f_gl;

% approx relative L^2-norm error vs the manufactured analyapprox)tical solution
xx = linspace(0,1,2*ncp);
approx = bspeval(p,[0 sol' 0],knot,xx);
exact = sin(2*pi*xx);

error = sqrt(sum((exact-approx).^2)/sum(exact.^2));

return