%clear b
clear

m=10; 
hh=1/m;
x=(0:hh:1)';     % assume uniform node spacing
for i=1:m,
  h(i)=x(i+1)-x(i);
end
h=h';

for i=1:m-1;
  b(i)=(quad8('xxf',x(i),x(i+1)) ...
              -x(i)*quad8('ff',x(i),x(i+1)))/h(i);
  b(i)=b(i)+(x(i+2)*quad8('ff',x(i+1),x(i+2)) ...
              - quad8('xxf',x(i+1),x(i+2)))/h(i+1);
end;
b=b';

% assume again uniform h
mdiag=diag(ones(m-1,1));
udiag=diag(ones(m-2,1),+1);
ldiag=diag(ones(m-2,1),-1);

lmat=-1/hh*ldiag+2/hh*mdiag-1/hh*udiag;

sol=[0; lmat\b; 0];

plot(x,sol,'d-'); hold on;
plot(0:0.01:1,u(0:0.01:1))
hold off;

% I'll skip the error calculation for now...
