function [y,x] = bvp_sd_dd(dgl, x0,xf,y0,yf,n);
% y = bvp_sd_dd(dgl, x0,xf,y0,yf,n);
% May 8, 2006
% Solution of 1 d 2 point bvp problem, symmetric differences, Dirichlet
% boundary conditions

h = (xf - x0)/n; h12 = h/2; h2 = h^2;
x =[x0+h:h:xf-h]';

[b,c,f] = feval(dgl,x);

A = diag(c*h2 -2) + diag(1+b(1:end-1)*h12,1) + diag(1-b(2:end)*h12,-1);

F = f*h2;
F(1) = F(1) - (1-b(1)*h12)*y0;
F(end) = F(end) - (1 + b(end)*h12)*yf;

y = A\F;