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

// *** 多項式(ラグランジュ)補間 ***
Pn = zeros(X);
for m = 1:length(X)
  for i = 1:(N + 1)
    Z(i) = 1;
    for j = 1:(N + 1)
      if i <> j
        Z(i) = Z(i) * (X(m) - Xn(j)) / (Xn(i) - Xn(j));
      end
    end
    Pn(m) = Pn(m) + Fn(i) * Z(i);
  end
end

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

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