function u=col_bvp2(n,x,f,p,q,am,bm,cm,ap,bp,cp)
%
% Solve -u_xx +p(x) u_x +q(x) u=f(x),  x in (-1,1)
% with
% am u(-1)+bm u'(-1)=cm, ap u(1)+bp u'(1)=cp,
%
% by using the spectral-collocation method
%
% Input: n, x(1:n+1): collocation points
%        am,bm,cm,ap,bp,cp: boundary conditions as described above
  %        f,p,q (1:n+1):  values of f, p, q at x(1:n+1)
% Output: u(1:n+1): the solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute the first and second derivative matrix

  dm=derivmatrix(n,2,x);

% form the big matrix

for j=1:n+1
    for i=1:n+1
	    a(i,j)=-dm(i,j,2)+p(i)*dm(i,j,1);
    end
    a(j,j)=a(j,j)+q(j);
end

% replace the first-row and last row by boundary conditions

    a(1,1)=am+bm*dm(1,1,1);
    a(n+1,n+1)=ap+bp*dm(n+1,n+1,1);
for i=2:n+1
     a(1,i)=bm*dm(1,i,1);
     a(n+1,i-1)=bp*dm(n+1,i-1,1);
end

% setup the righthand side
g(1)=cm;
g(n+1)=cp;
for i=2:n
	g(i)=f(i);
end
u=a\g';