clear;

// *** 近似する元の関数 ***
function y = f(x)
// y = 2 .* exp(x  - 1) - 1
 y = 2 ./ (1 + 9 .* x .^ 2) - 1
endfunction

// *** データ点 ***
N = 10;				// データ数
Xn = linspace(-1, 1, N + 1);	// xのデータ点
Fn = f(Xn);			// yのデータ点
X  = linspace(-1, 1, 100);	// プロット用のx

// *** 多項式(ラグランジュ)補間 ***
Wn = poly(0,"x");
Pn = 0;
for i = 1:length(Xn) do
  Wn = poly(Xn(find(Xn <> Xn(i))),"x");
  Pn = Pn + Fn(i) * Wn / horner(Wn,Xn(i));
end

PN = horner(Pn, X);

// *** グラフのプロット ***
subplot(2,1,1);
plot(X, f(X));
plot(Xn, Fn, 'sr');
plot(X,  PN, '--k');

// *** 誤差のプロット ***
subplot(2,1,2);
plot(X, PN - f(X));