n=20;
hh=rand(n-1,1)+0.01;
h=hh/sum(hh);
for i=1:n-1
   x(i+1)=x(i)+h(i)
end

A=zeros(18);
for i=1:18
    for j=1:18
        if i==j
            A(i,j)=1/h(i)+1/h(i+1);
        elseif j-i==1
            A(i,j)=-1/h(j);
        elseif i-j==1
            A(i,j)=-1/h(i);
        end
    end
end
A

for i=1:n-2
S(i)=sin(pi*x(i))*h(i)*(1/2)+sin(pi*x(i+1))*h(i)*(1/2)+sin(pi*x(i+1))*h(i+1)*(1/2)+sin(pi*x(i+2))*h(i+1)*(1/2);
end
S = 1/2 * pi^2*S;
u=inv(A)*S';
u = [0 u' 0];

plot(x,sin(pi*x), 'bo-',x,u,'r*-')