clear all;  

alf = 1/3; a=0; b=1; n=20; 
bl=1; br=1; 
h=(b-a)/n; x=a+h:h:b-h; x=x'; h1 =h*h; 
 
C=zeros(n-1,n-1); f = zeros(n-1,1);   
 
for i=2:n-2 
   beta_l = beta_h(x(i)-0.5*h-alf,bl,br,alf,h); 
   beta_r = beta_h(x(i)+0.5*h-alf,bl,br,alf,h); 
   C(i,i-1) = beta_l/h1;;C(i,i+1)=beta_r/h1; 
   C(i,i) = -(beta_l+beta_r)/h1; 
   f(i)=delta(x(i)-alf,h); 
end 
C(n-1,n-1)=-2*br/h1; C(n-1,n-2)=br/h1;C(1,1)=-2*bl/h1; C(1,2)=bl/h1; 
 
U=C\f; 
for i=1:n-1 
   soln(i)=uexact(x(i),bl,br,alf); e(i)=soln(i)-U(i); 
end 
figure(3);  plot(x,U,'o',x,soln);   
figure(4); plot(x,e)