NB. Find spline 2nd Derivative Matrix for either 1st 
NB. derivatives fixed at either or both endpoints or 
NB. free: xp0=[dy/dx] @ x_0,  xpN=[dy/dx] @ x_N 
NB. If one or both boundaries are free (y''=0),
NB. set derivative(s) xp0 and/or xpN = 'f'.
NB. Usage: (xp0;xpN) SpMDD x ---> B;C 
NB. Boxed output is a matrix B and vector C such that 
NB. (B mm y) +C --> D ,the spline fit 2nd derivatives
NB. for y(x), i.e., D = (xp0;x;xpN) splineD y 
NB. If both xp0 & xpN =0 or 'f', C will be all zeros.
NB. 
SpMDD=: 4 : 0
'xp0 xpN'=. x      NB. 1st derivatives on boundaries
k1=. 6% h1=. }: h=. DEL y 
k3=. 6% h3=. }. h
if. xp0 = 'f' do.  NB.========NB. first free &:
  if. xpN='f' do.  NB.--------NB. both free
    h2=. +: h3+h1
    h1=. 0,(}.h1)
    h3=. (}: h3),0
    AI=. trivert h1;h2;h3     NB. invert (N-2)x(N-2)
    AI=. 0,(|: 0,(|: AI),0),0 NB. border with zeros
    k1=. 0,k1,0
    k3=. 0,k3,0
    k2=. - (k1 + k3)
    B=. AI mm full_matrix k1;k2;k3
    B;(#y)#0 
  else.  NB.------------------NB. or just first free
    h2=. (+: h3+h1),(2*hN=.{:h)
    h1=. 0,(}.h1),hN
    h3=. h3,0
    AI=. trivert h1;h2;h3     NB. invert (N-1)x(N-1)
    AI=. 0,(|: 0,(|: AI))     NB. 0 border left & top 
    k2=. -(k1 + k3)
    k1=. 0,k1,(6%hN)
    k3=. 0,k3,0
    k2=. 0,k2,(_6%hN)
    B=. AI mm full_matrix k1;k2;k3
    B; AI mm ((#h)#0),(6*xpN)  
  end.
else.    NB.==================NB. or first derivative set
  if. xpN='f' do. NB.---------NB. & last free
    h2=. (2*h0=.{.h),(+: h3+h1)
    h1=. 0,h1
    h3=. h0,h3
    AI=. trivert h1;h2;h3     NB. invert (N-1)x(N-1)
    AI=. (|: (|: AI),0),0     NB. 0 border right & bottom 
    k2=. -(k1+k3)
    k1=. 0,k1,0
    k3=. (6%h0),k3,0
    k2=. (_6%h0),k2,0
    B=. AI mm full_matrix k1;k2;k3
    B; AI mm (_6*xp0),((#h)#0) 
  else.  NB. -----------------NB. or both derivatives set
    h2=. (2*h0=.{.h),(+: h3+h1),(2*hN=.{:h)
    h1=. 0,h1,hN
    h3=. h0,h3,0
    AI=. trivert h1;h2;h3     NB. invert NxN 
    k2=. -(k1 + k3)
    k1=. 0,k1,(6%hN)
    k3=. (6%h0),k3,0
    k2=. (_6%h0),k2,(_6%hN)
    B=. AI mm full_matrix k1;k2;k3
    B; AI mm (_6*xp0),((<:#h)#0),(6*xpN)
  end.
end.
) 

NB. Invert a tridiagonal matrix by direct substitution.
NB. Rows are  a_i*x_(i-1), b_i*x_i, c_i*x_(i+1)  
NB. a,b,c have same length, i.e., a[0]=0 and c[n]=0. 
NB. Usage: trivert a;b;c --> inverse(full nxn matrix) 
NB. 
trivert=: 3 : 0
'a b c'=. y
n1=. <: n=. <: nn=. #b
x=. %{.b
inv=. x (<0 0)}(=@i. nn)  
v=. (x*{.c) 0}v=. nn#0
j=. 0
while. n1>:j=.>:j do.
  b1=. j{b
  a1=. j{a
  v0=. (<:j){v
  x=. %b1-a1*v0
  c1=. j{c
  v=. (x*c1)j}v
  inv1=. j{inv
  inv0=. (<:j){inv
  inv=. (x*inv1-a1*inv0)j}inv
end.
x=. %(n{b)-(n{a)*(n1{v)
inv=. (x*(n{inv)-(n{a)*(n1{inv)) n}inv
j=. n
while. 0<:j=.<:j do.
  inv1=. j{inv
  v1=. j{v
  inv2=. (>:j){inv
  inv=. (inv1-v1*inv2)j}inv
end.
)

NB. Expand the three diagonal bands of a tridiagonal 
NB. matrix into the full matrix. 
NB. Rows are  a_i, b_i, c_i
NB. a,b,c have same length, i.e., a[0]=0 and c[n]=0.
NB. Usage: full_matrix a;b;c --> (full nxn matrix)
NB.
full_matrix=: 3 : 0
'a b c'=. y
abc=. b* I=. =@i.(#b)   NB. I is the identity matrix
abc=. abc+ 1|."1 (a* I)
abc+  _1|."1 (c* I)
)