load '/your_path/MTs.ijs' 
NB. (matrix transforms for slab geometry)
load '/your_path/s_and_p.ijs' 
load '/your_path/IQMs.ijs' 
load '/your_path/IQs.ijs' 

NB. define tau grid: 1st point (after 0) is 0.01, last
NB. point is 2, logarithmic spacing, 8 points per decade:
tau=: tau_grid2 0.01 8 1.0
NB. 35 points: 0 0.01 0.0133352 0.0177828... 1.98666 1.99 2

L0=. 1 MTs tau       NB. get Lam, M & N transforms
M0=. 3 MTs tau
N0=. 5 MTs tau

NB. Example of solution for lambda=0.5 and parabolic B(tau):
B=. 1-(tau-1)^2                    NB. parabolic source term
lam=. (#B)#0.5                     NB. constant lambda
's p'=: (L0;M0;N0) s_and_p (lam;B) NB. solve for s & p
load 'plot'                        NB. load J plotting package
'pensize 2' plot tau;>B;s;10*p     NB. plot B, s and 10 p vs. tau
NB. the above plot will get overwritten quickly 
NB. define a grid of mu (= cos(theta)) values:
mu=: int01 40
I0=. IQMs tau;mu              NB. get matrix for I(0,mu)
'I Q pol'=: (I0;mu) IQs s;p   NB. evaluate emergent radiation
NB. plot I(mu) and -10 times the polarization as function of mu
'pensize 2' plot mu;>I;_10*pol

NB. Construct logarithmic grid over the interval [0,2T],
NB. symmetric about central palne at tau=T.
NB. tau1 <- first optical depth point, e.g. tau1=: 0.01
NB. npd <- number of points/decade, e.g. npd=: 8
NB. T <- Half-depth of slab, e.g. T=: 1
tau_grid2=: 3 : 0
'tau1 npd T'=: y
t=. T,~(i. +/t<T){ t=: 0,tau1*10^(npd %~ i. 500)  NB. tau points
t, }. |. (2*T)-t
)