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*-')