--- a/scripts/polynomial/spline.m	Thu Apr 19 17:44:46 2012 -0400
+++ b/scripts/polynomial/spline.m	Wed Apr 11 12:54:41 2012 +0200
@@ -131,31 +131,12 @@
       b = b(1:n-1,:);
       a = a(1:n-1,:);
     else
-      if (n == 3)
-        dg = 1.5 * h(1) - 0.5 * h(2);
-        c(2:n-1,:) = 1/dg(1);
-      else
-        dg = 2 * (h(1:n-2) .+ h(2:n-1));
-        dg(1) = dg(1) - 0.5 * h(1);
-        dg(n-2) = dg(n-2) - 0.5 * h(n-1);
-
-        e = h(2:n-2);
-
-        g = 3 * diff (a(2:n,:)) ./ h(2:n-1,idx) ...
-          - 3 * diff (a(1:n-1,:)) ./ h(1:n-2,idx);
-        g(1,:) = 3 * (a(3,:) - a(2,:)) / h(2) ...
-          - 3 / 2 * (3 * (a(2,:) - a(1,:)) / h(1) - dfs);
-        g(n-2,:) = 3 / 2 * (3 * (a(n,:) - a(n-1,:)) / h(n-1) - dfe) ...
-          - 3 * (a(n-1,:) - a(n-2,:)) / h(n-2);
-
-        c(2:n-1,:) = spdiags ([[e(:); 0], dg, [0; e(:)]],
-                              [-1, 0, 1], n-2, n-2) \ g;
-      endif
-
-      c(1,:) = (3 / h(1) * (a(2,:) - a(1,:)) - 3 * dfs
-             - c(2,:) * h(1)) / (2 * h(1));
-      c(n,:) = - (3 / h(n-1) * (a(n,:) - a(n-1,:)) - 3 * dfe
-             + c(n-1,:) * h(n-1)) / (2 * h(n-1));
+      g(1,:) = (a(2,:) - a(1,:)) / h(1) - dfs;
+      g(2:n-1,:) = (a(3:n,:) - a(2:n-1,:)) ./ h(2:n-1) - ...
+                   (a(2:n-1,:) - a(1:n-2,:)) ./ h(1:n-2);
+      g(n,:) = dfe - (a(n,:) - a(n-1,:)) / h(n-1);
+      c = spdiags([[h/6;0],[h(1)/3;(h(1:n-2)+h(2:n-1))/3;h(n-1)/3],[0;h/6]],...
+                  [-1,0,1],n,n) \ (g / 2);
       b(1:n-1,:) = diff (a) ./ h(1:n-1, idx) ...
         - h(1:n-1,idx) / 3 .* (c(2:n,:) + 2 * c(1:n-1,:));
       d = diff (c) ./ (3 * h(1:n-1, idx));
@@ -297,11 +278,23 @@
 %! y = [1,2,3,4];
 %! pp = spline (x,y);
 %! [x,P] = unmkpp (pp);
-%! assert (norm (P-[3,-3,1,2]), 0, abserr);
+%! assert (P, [3,-3,1,2], abserr);
 %!test
 %! x = [2,1];
 %! y = [1,2,3,4];
 %! pp = spline (x,y);
 %! [x,P] = unmkpp (pp);
-%! assert (norm (P-[7,-9,1,3]), 0, abserr);
+%! assert (P, [7,-9,1,3], abserr);
+%!test
+%! x = [0,1,2];
+%! y = [0,0,1,0,0];
+%! pp = spline (x,y);
+%! [x,P] = unmkpp (pp);
+%! assert (P, [-2,3,0,0;2,-3,0,1], abserr);
+%!test
+%! x = [0,1,2,3];
+%! y = [0,0,1,1,0,0];
+%! pp = spline (x,y);
+%! [x,P] = unmkpp (pp);
+%! assert (P, [-1,2,0,0;0,-1,1,1;1,-1,-1,1], abserr);