+
+ alpha = (double *)malloc(sizeof(double) * (n-2));
+ for (i=0; i<n-2; i++) {
+ alpha[i] = 2 * (h[i] + h[i+1]);
+ }
+
+ /* B[] is the vector on the right hand side of the equation */
+ B = (double *)malloc(sizeof(double) * (n-2));
+ for (i=0; i<n-2; i++) {
+ B[i] = 6 * ((f[i+2] - f[i+1])/h[i+1] - (f[i+1] - f[i])/h[i]);
+ }
+
+ /* Now solve the n-2 by n-2 system */
+ m = (double *)malloc(sizeof(double) * (n-2));
+ for (i=1; i<=n-3; i++) {
+ /*
+ d0 = alpha 0
+ a0 = h1
+ c0 = h1
+
+ di = di - (ai-1 / di-1) * ci-1
+ bi = bi - (ai-1 / di-1) * bi-1
+ ;
+ */
+ alpha[i] = alpha[i] - (h[i]/alpha[i-1]) * h[i];
+ B[i] = B[i] - (h[i]/alpha[i-1]) * B[i-1];
+ }
+ /* xn-3 = bn-3 / dn-3; */
+ m[n-3] = B[n-3]/alpha[n-3];
+ for (i=n-4; i>=0; i--) {
+ m[i] = (B[i]-h[i+1]*m[i+1])/alpha[i];
+ }
+
+ for (i=0; i<orig_n-1; i++) {
+ double mi, mi1;
+ mi = (i==(n-2)) ? 0 : m[i];
+ mi1 = (i==0) ? 0 : m[i-1];
+
+ p[i].a = f[i+1];
+ p[i].b = (f[i+1] - f[i]) / h[i] + h[i] * (2*mi + mi1) / 6;
+ p[i].c = mi/2;
+ p[i].d = (mi-mi1)/(6*h[i]);
+ }
+
+ free(alpha);
+ free(B);
+ free(h);
+ free(m);