NB. Matrix transforms for symmetric slab
NB. idx=1,2,3,4,5 -> Lambda,Phi,M,Phi^4,N transforms
NB. idx MT tau_grid -> (linear);(spline) matrix operator
NB. The quantities ending in "f" refer to the far side.
load '/your_path/EIn.ijs'
load '/your_path/SpMDD.ijs'

MTss=: 4 : 0
pm=. }:"(1) _1+ +: t<:/t=. y
tf=. (+: {: t)-t  NB. far side:  2*tau(N)-tau(i)
'I0 If0 I1 If1 I2 If2 I3 If3'=. x IIss t
t3=. t* t2=. t*t
tf3=. tf* tf2=. tf*tf
a2=. *: a=. }: t
b2=. *: b=. }. t
h2=. *: h=. b-a
h6=. 6%~ hi=. % h
X=. hi*"1 ((-t-/b)*I0) -pm*I1
Xf=. hi*"1 ((-tf-/b)*If0) +If1
Y=. -hi*"1 ((-t-/a)*I0) -pm*I1
Yf=. -hi*"1 ((-tf-/a)*If0) +If1
lin=. (X,"(1) 0)+(0,"(1) Y)
lin=. lin+ (Xf,"(1) 0)+(0,"(1) Yf)
NB. 1's, length ur= #rows, uc= #columns:
uc=. }. ur=. (#t)#1 
u1=. ur*/ b*(b2-h2)
u2=. -t*/((3*b2)-h2)
us2=. -tf*/((3*b2)-h2)
u3=. ((3*t2)*/b)-(t3*/uc)
us3=. ((3*tf2)*/b)-(tf3*/uc)
U=. I0* u1+u2+u3
US=. If0* u1+us2+us3
u1=. -ur*/ (3*b2)-h2
u2=. (6*t*/b)-(3*t2)*/uc
us2=. (6*tf*/b)-(3*tf2)*/uc
U=. U+ (pm*I1*(u1+u2))
US=. US -If1*(u1+us2)
U=. U+ (I2*3*(-t-/b))
US=. US+ If2*3*(-tf-/b)
U=. h6*"1 U -pm*I3
U=. U+ h6*"1 US+ If3
v1=. ur*/ a*(a2-h2)
v2=. -t*/((3*a2)-h2)
vs2=. -tf*/((3*a2)-h2)
v3=. (3*t2*/a)-(t3*/uc)
vs3=. (3*tf2*/a)-(tf3*/uc)
V=. -I0* v1+v2+v3 
VS=. -If0* v1+vs2+vs3 
v1=. ur*/ (3*a2)-h2
v2=. ((3*t2)*/uc)-6*t*/a
vs2=. ((3*tf2)*/uc)-6*tf*/a
V=. V+ pm*I1*(v1+v2)
VS=. VS- If1*(v1+vs2)
V=. V+ I2*3*(t-/a) 
VS=. VS+ If2*3*(tf-/a)
V=. h6*"1 V+ pm*I3
V=. V+ h6*"1 VS- If3
C=. >0{ ('f';0) SpMDD t
spl=. (U mm }:C)+(V mm }.C)
lin+spl
)

NB. make table of I-integrals for MTss  
NB. symmetric slab: tf=tau to far side       
IIss=: 4 : 0
idx=. x
a=. | y -/ (}:y) NB. |tau(k)-tau(i)|
b=. | y -/ (}.y) NB. |tau(k+1)-tau(i)|
pm=. }:"(1) _1+ +: y<:/y
if. -.(idx=2)+:(idx=4) do. pm=. |pm end.
tf=. (+: {: y)-y NB. 2*tau(N)-tau(i)
c=. tf -/ (}:y)  NB. tf(k)-tau(i)
d=. tf -/ (}.y)  NB. tf(k+1)-tau(i)
if. idx=1 do.
NB. the kernels for Lambda transform
K0=. 3 : '_2%~ 2&EIn y'
K1=. 3 : '2%~ (3&EIn y)-(^-y)'
K2=. 3 : '-(4&EIn y)+ 2%~y*(^-y)'
K3=. 3 : '(3*5&EIn y)- 2%~(^-y)*(3+y+y*y)'
end.
if. idx=2 do.
NB. the kernels for Phi-transform
K0=. 3 : '_2*3&EIn y'
K1=. 3 : '+:(2*4&EIn y)- ^-y'
K2=. 3 : '(_12*5&EIn y)+2*(^-y)*(1-y)'
K3=. 3 : '(48*6&EIn y)-2*(^-y)*(6+y*y)'
end.
if. idx=3 do.
NB. the kernels for M-transform
K0=. 3 : '(_2%~ 2&EIn y)+(3%2)*4&EIn y'
K1=. 3 : '-:(3&EIn y)+(_9*5&EIn y)+2*(^-y)'
K2=. 3 : '(-(4&EIn y))+(18*6&EIn y)-(^-y)*(3-y)'
K3=. 3 : '(3*5&EIn y)+(_90*7&EIn y)+(^-y)*(15+(_2*y)+y*y)'
end.
if. idx=4 do.
NB. the kernels for Phi^(4)
K0=. 3 : '((_2%3)*3&EIn y)+(2*5&EIn y)'
K1=. 3 : '((4%3)*4&EIn y)+(_8*6&EIn y)+(4%3)*(^-y)'
K2=. 3 : '(_4*5&EIn y)+(40*7&EIn y)-(4%3)*(^-y)*(4-y)'
K3=. 3 : '(16*6&EIn y)+(_240*8&EIn y)+(4%3)*(^-y)*(24+(_3*y)+y*y)'
end.
if. idx=5 do.
NB. the kernels for N_tau:
K0=. 3 : '((_5%3)*2&EIn y)+(4*4&EIn y)-(3*6&EIn y)'
K1=. 3 : '((5%3)*3&EIn y)+(_12*5&EIn y)+(15*7&EIn y)-(2%3)*(^-y)'
K2=. 3 : '((_10%3)*4&EIn y)+(48*6&EIn y)+(_90*8&EIn y)+(2%3)*(^-y)*(6-y)'
K3=. 3 : '(10*5&EIn y)+(_240*7&EIn y)+(630*9&EIn y)-(2%3)*(^-y)*(63+(_5*y)+y*y)'
end.
I0=. pm*(K0 b) -K0 a
If0=. (K0 c) -K0 d
I1=. pm*(K1 b) -K1 a
If1=. (K1 c) -K1 d
I2=. pm*(K2 b) -K2 a
If2=. (K2 c) -K2 d
I3=. pm*(K3 b) -K3 a
If3=. (K3 c) -K3 d
I0;If0;I1;If1;I2;If2;I3;If3
)