clear;clc;
n=20;
h=rand(n-1,1)+0.01;
h=h/sum(h);
for i=1:n-1
x(1)=0;
x(i+1)=x(i)+h(i);
end
A=zeros(n-2,n-2);
for i=1:n-2
for j=1:n-2
if (i==j)
A(i,j)=1/h(i)+1/h(i+1);
elseif (i-j==-1)
A(i,j)=-1/h(j);
elseif (i-j==1)
A(i,j)=-1/h(i);
end
end
end
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
ss = 1/2 * pi^2*s;
u=inv(A)*ss';
u = [0 u' 0];
plot(x,sin(pi*x), 'ro-')
hold on
plot(x,u,'k*-')