function u=col_bih1d(n,x,f,p,q,am,bm,ap,bp)
%
% Solve u_xxxx -p(x) u_xx +q(x) u=f(x),  x in (-1,1)
% with
% u(-1)=am, u'(-1)=bm,  u(1)=ap, u'(1)=bp
%
% by using the spectral-collocation method
%
% Input: n, x(1:n+1): collocation points
%        am,bm,ap,bp: boundary conditions as described above
%        f,p,q: external functions
% Output: u(1:n+1): the solution
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

% compute the first and second derivative matrix

  dm=derivmatrix(n,4,x);

% form the big matrix

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

% replace the first two rows and last two rows by boundary conditions

    a(1,1)=1;
a(n+1,n+1)=1;
for i=2:n+1
a(1,i)=0;
a(n+1,i-1)=0;
end
for i=1:n+1
a(2,i)=dm(1,i,1);
a(n,i)=dm(n+1,i,1);
end
% setup the righthand side
g(1)=am;
g(2)=bm;
g(n)=bp;
g(n+1)=ap;
for i=3:n-1
	g(i)=feval(f,x(i));
end
u=a\g';
cond(a)