clear all; close all; clf
N=636;
x=linspace(0,1,N);
h=rand(N-1,1)+636;
h=h/sum(h);
for i=1:N-1
x(i+1)=x(i)+h(i);
end
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 || i==j-1
A(i,j)=-1/h(max(i,j));
end
end
end
for i=1:N-2
s(i) = h(i)*sin(pi*x(i))/2 + h(i)*sin(pi*x(i+1))/2 + ...
h(i+1)*sin(pi*x(i))/2 + h(i+1)*sin(pi*x(i+1))/2;
end
s=pi^2*s/2;
u=inv(A)*s';
u=[0 u' 0];
for i=1:N-2
y(i)=hat(x(i+1),x(i),x(i+1),x(i+2));
end
y = [0 y 0];
for i=1:N
f(i)=(u(i)*y(i));
end
plot(x,sin(pi*x),'r-',x,f,'b--')
grid on