function [x,u] =ode1d(a,b,aa,bb,n)
% solve -y"(x) + p(x) y'(x) + q(x)y(x) = r(x)
% on a grid of n +2 points from a to b
% Dirichlet b.c.s : y(a) = aa, y(b) = bb
h = (b-a)/(n+1);   % grid spacing
x = [a + h : h : b-h];  % n samples
A = x; B = x; C = x; S = x;
for i = 1:n
% create the "full" diagonals
A(i) = -0.5*(1 + h/2*p(x(i)));
B(i) = (1 + h*h/2*q(x(i)));
C(i) = -0.5*(1 - h/2*p(x(i)));
% the right hand side
S(i) = h*h/2*r(x(i));
end
% and the boundary conditions
S(1) = S(1) + 0.5*(1 + h/2*p(a))*aa;
S(n) = S(n) + 0.5*(1 - h/2*p(b))*bb;

% make a sparse matrix
B = [A',B',C'];
d = [-1,0,1];
%S
K = spdiags(B,d,n,n);
%full(K);
u = K\S';