[跪求][帮忙看看]写了个C# 三次样条插值的类库-三转角法(自然边界)

各位帮忙看看, 第一次做这种数值类计算,完全按照公式逐步求解, 可是当我发现一个严重问题
就是加入原函数为x=2y, 当在 x 上10以后的区间插值就插值出来的曲线就开始剧烈波动了.....
各位大侠帮忙看看其中到底哪里出了问题:
using System;
using System.Collections.Generic;
using System.Text;namespace triSpline
{
    public class triSpline
    {
        /// <summary>
        /// 三次插值(三转角法)
        /// </summary>
        /// <param name="x">节点x值</param>
        /// <param name="y">节点y值</param>
        /// <param name="n">节点数量</param>
        /// <param name="t_x">插值节点x列表</param>
        /// <param name="m">插值节点数量</param>
        /// <returns>插值节点y列表</returns>
        public static double[] GettriSpline_y(double[] x, double[] y, int n, double[] t_x, int m)
        {
            double[] t_y = new double[m]; //保存插值结果y的列表            #region 数据解说
            //原输入节点为x[0]~x[n]: 
            // 2*m[0] + m[1] =d[0]
            // λ[j]*m[j-1] + 2*m[j] + μ[j]*m[j+1] = d[j]
            // m[n-1] + 2*m[n] = d[n]
            //
            //当前输入为x[0]~x[n-1]，则应为:
            // 2*m[0] + m[1] =d[0]
            // λ[j]*m[j-1] + 2*m[j] + μ[j]*m[j+1] = d[j]
            // m[n-2] + 2*m[n-1] = d[n-1]            //设定解LU时用的系数列表,即三对角线的矩阵(对应线性方程组系数矩阵)中 a=lamda,b=2,c=ita
            #endregion            double[] b = new double[n];     //(即b-2)保存三对角线中间数值列表
            double[] a = new double[n]; //(即a-λ lamda )由于最后用来保存一阶导数值 需要n个, 所以最后一个a[n-1]在前期空出来
            double[] c = new double[n - 1]; //(即c-μ ita)最后一行a[n-1]没有c, c在对角线上侧
            double[] d = new double[n];     //(即f or g or 方程中的d) 应该有n个,第一个和最后一个是特殊的, 中间n-2个具有通项
            double[] hx = new double[n - 1]; //保存x[j+1]-x[j] 
            double[] hy = new double[n - 1];   //保存y[j+1]-y[j] 用来计算f[x0，x1]等
            double[] f = new double[n - 1];    //f[x0，x1]
            #region 原循环, 从0开始到n-2 共n-1次
            ////先解决a b c 表达式, 用于追赶法求解 LU 的 Alpha,Beta 
            //for (int j = 0; j < n - 1; ++j) //0~n-2循环  最后一次未赋值,比输入数值n少1个,共n-1次循环,
            //{
            //    b[j] = 2; //中对角线为2值不变 n             //    hx[j + 1] = x[j + 1] - x[j]; //n-1个 比输入数值少1个
            //    hy[j + 1] = x[j + 1] - x[j];
            //    f[j+1] = hy[j + 1] / hx[j + 1];            //    if (j > 0) //第一个元素为特殊值,由此判定跳过
            //    {
            //        a[j - 1] = hx[j + 1] / (hx[j + 1] + hx[j]); // 由于a为数组必须向前移一位从a[0]开始赋值,但总数只有n-1个
            //        c[j] = 1 - a[j - 1];
            //        d[j] = 3 * (a[j - 1] * f[j] + c[j] * f[j+1]);//可带入第一个和最后一个确认
            //    }
            //}
            #endregion            //补充第一次循环
            hx[0] = x[1] - x[0]; //末元素[n-2] 共n-1个元素
            hy[0] = y[1] - y[0];  //末元素[n-2] 共n-1个元素
            f[0] = hy[0] / hx[0]; //末元素[n-2] 共n-1个元素
            //补充矩阵第一行(第一行无a-λ参数)
            b[0] = 2;
            c[0] = 1;    //第一行的c值
            d[0] = 3 * f[0]; //因为自然条件 -h/2*f"=0            for (int j = 1; j < n - 1; ++j) //1~n-2循环 去掉最后一行和第一行,由此去掉if判断
            {
                b[j] = 2; //中对角线为2值不变 n                 hx[j] = x[j + 1] - x[j]; //n-1个 比输入数值少1个
                hy[j] = y[j + 1] - y[j];
                f[j] = hy[j] / hx[j];                a[j - 1] = hx[j] / (hx[j] + hx[j - 1]); // 由于a为数组必须向前移一位从a[0]开始赋值,但总数只有n-1个
                c[j] = 1 - a[j - 1];
                d[j] = 3 * (a[j - 1] * f[j - 1] + c[j] * f[j]);//可带入第一个和最后一个确认   
            }            //补充矩阵最后一行特殊值
            b[n - 1] = 2;  //最后一个值,以为上述循环未给最后一个赋值
            a[n - 2] = 1;  //最后一行的a值,由于a为数组,第一个a[0]不能为空,所以最后n-1向前移动一位变为n-2
            d[n - 1] = 3 * f[n - 2]; //因为自然条件 -h/2*f"=0            //分解计算公式 A=L*U
            for (int i = 0; i < n - 1; ++i)
            {
                c[i] /= b[i]; //将Beta 保存在c[]中 0~n-2
                b[i + 1]  -= a[i] * c[i];//将Alpha 保存在b[]中 1~n-1, b[0]与Alpha[0]相同,所以不做变动.
            }            //解方程组Ly=f, 将求得的y存于b[]中, 共0~n-1个值 即n个值
            b[0] = y[0] / b[0];
            for (int i = 1; i < n; ++i)
            {
                b[i] = (y[i] - a[i - 1] * b[i - 1]) / b[i];
            }            //解Ux=y,将求得的x存于a[]中,最后一个x 保存在a[n-1]中, 剩下循环0~n-2, 所以此时a[] 共n个值
            a[n - 1] = b[n - 1];
            for (int i = n - 2; i >= 0; --i)
            {
                a[i] = b[i] - c[i] * a[i + 1];
            }
            for (int j = 0; j < m; ++j) //可外插大于边缘的数值(不推荐使用)
            {
                int i;
                if (t_x[j] >= x[n - 1])
                    i = n - 2;
                else
                {
                    i = 0;
                    while (t_x[j] > x[i + 1])
                        ++i;
                }
                    t_y[j] = y[i] * Math.Pow(x[i + 1] - t_x[j], 2) * (hx[i] + 2 * (t_x[j] - x[i])) / Math.Pow(hx[i], 3)
                            + y[i + 1] * Math.Pow(x[i] - t_x[j], 2) * (hx[i] + 2 * (x[i + 1] - t_x[j])) / Math.Pow(hx[i], 3)
                            + a[i] * Math.Pow(x[i + 1] - t_x[j], 2) * (t_x[j] - x[i]) / Math.Pow(hx[i], 2)
                            + a[i + 1] * Math.Pow(x[i] - t_x[j], 2) * (t_x[j] - x[i + 1]) / Math.Pow(hx[i], 2);
                //}
            }
            return t_y;
        }
    }
}