lib/interpolation/one_dimensional.rb
def compute_second_derivatives_for y
      y_sd = Array.new(@size)

      n      = y_sd.size
      u      = Array.new(n-1)
      yp1    = @opts[:yp1] # first derivative of the 0th point as specified by the user
      ypn    = @opts[:ypn] # first derivative of the nth point as specified by the user
      qn, un = nil, nil

      if yp1 > 0.99E30
        y_sd[0], u[0] = 0.0, 0.0
      else 
        y_sd[0] = -0.5
        u[0] = (3.0 / (@x[1] - @x[0])) * ((y[1] - y[0]) / (@x[1] - @x[0]) - yp1)
      end
      
      1.upto(n-2) do |i| # decomposition loop for tridiagonal algorithm
        sig      = ( @x[i] - @x[i-1] ) / ( @x[i+1] - @x[i-1] )
        p        = sig * y_sd[i-1] + 2
        y_sd[i] = ( sig - 1) / p 
        u[i]     = (( y[i+1] - y[i]) / (@x[i+1] - @x[i])) - ((y[i] - y[i-1]) / (@x[i] - @x[i-1]))
        u[i]     = ( 6 * u[i] / ( @x[i+1] - @x[i-1] ) - sig * u[i-1] ) / p;
      end

      if ypn > 0.99E30
        qn, un = 0.0, 0.0
      else
        qn = 0.5
        un = (3.0 / ( @x[n-1] - @x[n-2] )) * ( ypn - ( y[n-1] - y[n-2] ) / ( @x[n-1] - @x[n-2] ))
      end
      y_sd[n-1] = ( un - qn * u[n-2] ) / ( qn * y_sd[n-2] + 1.0 )

      (n-2).downto(0) do |k|
        y_sd[k] = y_sd[k] * y_sd[k+1] + u[k]
      end

      y_sd
    end
def compute_second_derivatives_for y
      y_sd = Array.new(@size)

      n      = y_sd.size
      u      = Array.new(n-1)
      yp1    = @opts[:yp1] # first derivative of the 0th point as specified by the user
      ypn    = @opts[:ypn] # first derivative of the nth point as specified by the user
      qn, un = nil, nil

      if yp1 > 0.99E30
        y_sd[0], u[0] = 0.0, 0.0
      else 
        y_sd[0] = -0.5
        u[0] = (3.0 / (@x[1] - @x[0])) * ((y[1] - y[0]) / (@x[1] - @x[0]) - yp1)
      end
      
      1.upto(n-2) do |i| # decomposition loop for tridiagonal algorithm
        sig      = ( @x[i] - @x[i-1] ) / ( @x[i+1] - @x[i-1] )
        p        = sig * y_sd[i-1] + 2
        y_sd[i] = ( sig - 1) / p 
        u[i]     = (( y[i+1] - y[i]) / (@x[i+1] - @x[i])) - ((y[i] - y[i-1]) / (@x[i] - @x[i-1]))
        u[i]     = ( 6 * u[i] / ( @x[i+1] - @x[i-1] ) - sig * u[i-1] ) / p;
      end

      if ypn > 0.99E30
        qn, un = 0.0, 0.0
      else
        qn = 0.5
        un = (3.0 / ( @x[n-1] - @x[n-2] )) * ( ypn - ( y[n-1] - y[n-2] ) / ( @x[n-1] - @x[n-2] ))
      end
      y_sd[n-1] = ( un - qn * u[n-2] ) / ( qn * y_sd[n-2] + 1.0 )

      (n-2).downto(0) do |k|
        y_sd[k] = y_sd[k] * y_sd[k+1] + u[k]
      end

      y_sd
    end
References: Numerical Recipes Edition 3. Chapter 3.3