% function y = linfd(a,b,n,alpha,beta,p,q,r)
%
% Solves the linear BVP
%  y''(x) = p(x) y' + q(x) y(x) + r(x)
%  y(a)=alpha;  y(b)=beta
%
%
% Inputs
%  a,b    time interval
%
%  n      number of time subintervals
%
%  alpha  = y(a)
%
%  beta   = y(b)
%
%  p      function handle for p(x) (assumed vectorized)
%
%  q      function handle for q(x) (assumed vectorized)
%
%  r      function handle for r(x) (assumed vectorized)

function y = linfd(a,b,n,alpha,beta,p,q,r)
% setup grid
h=(b-a)/(n+1);  z = a + (1:n)'*h;
% evaluate at grid points
ps = p(z); qs = q(z); rs=r(z);
% setup matrix of linear system
L = (1/h^2)*spdiags(TO BE COMPLETED);
D = (1/2/h)*spdiags(TO BE COMPLETED);
Q = spdiags(TO BE COMPLETED);
A = L-D-Q;
% setup RHS of linear system
b = TO BE COMPLETED 

% assemble solution
y = TO BE COMPLETED
