function M = spline_coef(x, y)

  [~,n] = size(x);

  eq = zeros(n-2);

  for i=2:n-1
      h1 = x(i+1) - x(i);
      h0 = x(i) - x(i-1);
      lambda = h1/(h0+h1);
      mu = 1 - lambda;
      eq(i, i-1:i+1) = [mu 2 lambda];
      d(i) = 6 * ( (y(i+1) - y(i)) / h1 - (y(i) - y(i-1)) / h0 ) / (h0 + h1);
  end

  d(n) = 0;
  eq(1,1) = 2;
  eq(n,n) = 2;

  M = eq\d';

end