clear all;  

alf = 1/3; a=0; b=1; n=20; 
bl=1;  br=1;   bjump=br-bl; 
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=1:n-2 
   if  x(i) <=alf & alf< x(i+1) 
      j = i; 
   end 
end 
 
for i=2:n-2 
   tmp = my_beta(x(i),bl,br,alf); 
   C(i,i-1) = tmp/h1;C(i,i) = -2*tmp/h1;C(i,i+1)=tmp/h1; 
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; 
 
% At x=x(j) 
dj =h1+bjump*(x(j-1)-alf)*(x(j)-alf )/(2*bl); 
C(j,j-1) = (bl-bjump*(x(j)-alf )/h )/dj;C(j,j+1)=br/dj; 
C(j,j) = (-2*bl+bjump*(x(j-1)-alf )/h )/dj; 
f(j) = C(j,j+1)*(x(j+1)-alf)/br; 
 
dj =h1 - bjump*(x(j+1)-alf)*(x(j+2)-alf )/(2*br); 
C(j+1,j+2) = (br-bjump*(x(j+1)-alf )/h )/dj;C(j+1,j)=bl/dj; 
C(j+1,j+1) = (-2*br+bjump*(x(j+2)-alf )/h )/dj; 
f(j+1) = -C(j+1,j)*(x(j)-alf )/bl; 
 
U=C\f; 
for i=1:n-1 
   soln(i)=uexact(x(i),bl,br,alf); e(i)=soln(i)-U(i); 
end 
figure(1);  plot(x,U,'o',x,soln);   
figure(2); plot(x,e)