function [u,x]=fdBVP1d(a,b,c,f,u0,u1,M)

% function [u,x]=fdBVP1d(a,b,c,f,u0,u1,M)
% 
% FD for 1d BVP (notations in Larsson-Thomee sect 4.1)
%
% Input:
%     a,b,c   : coefficients of differential operator
%     f       : source term (rhs)
%     u0,u1   : (Dirichlet) boundary values
%     M       : number of intervals
% Output:
%     u       : solution (includes u(0) and u(1))
%     x       : mesh

h=1/M;
x=linspace(h,1-h,M-1);      % alt: x=linspace(0,1,M+1); x([1 end])=[];

% Matrix

a=a(x);       % continuous function substituted by grid function
b=b(x);
c=c(x);
d=2*a+h^2*c;                      % diagonal
u=-a(1:end-1)+(h/2)*b(1:end-1);   % upper diagonal
l=-a(2:end)-(h/2)*b(2:end);       % lower diagonal
i=[1:M-1   1:M-2   2:M-1];
j=[1:M-1   2:M-1   1:M-2];
matrix=sparse(i,j,[d u l]);             % (M-1) x (M-1)

% Right-hand side

rhs=h^2*f(x');
rhs(1)=rhs(1)+(a(1)+h/2*b(1))*u0;
rhs(end)=rhs(end)+(a(end)-h/2*b(end))*u1;

% Solution and post-processing

u=matrix\rhs;
u=[u0; u; u1];
x=[0;x';1];     

return