How to run an example of the Matrix Method for solution of radiative
transfer in a stellar atmosphere with Rayleigh (dipole) scattering.
This code is written in the "J" language (a decendent of APL). You can get
it at this URL: https://code.jsoftware.com/wiki/System/Installation/J9.7/Zips
(The latest beta.) It is free and cross-platform. 
First, download and install the J language and start J. To run this example, 
enter the commands given below. Please note that NB. (note bene) and any 
following text are just comments -- they do not need to be entered.

load 'Matrix_Method'      NB. load the J routines.
taus=. tau_grid ''        NB. define an 80 element tau grid (mostly log)
Set_LMN''                 NB. create matrix operators \Lamda, M and N
B=. 0.75*((2%3) + taus)   NB. B(tau): Eddington approx. to grey atmos.
llam=: {{0.01+0.99*y%(1+y)}} NB. arbitrary function for lam(tau)
lam=. llam taus           NB. set lam(tau) (scattering dominant near surface)
's p'=. (L0;M0;N0) s_and_p (lam ;B)  NB. solve for source vectors s & p
mu=. 0.05*i. 21           NB. define a grid of 21 values 0<= mu <=1
E=. IQMsi taus; mu        NB. create matrix operator E for emergent rays
'I Q P'=. (E;mu) IQs s;p  NB. evaluate Stokes I & Q, add P=Q/I polarization
load'plot'                NB. load J plot packkage
ps=. 'pensize 3'          NB. make lines thicker
ps plot mu; _100*P        NB. plot percent polarization vs. mu (mu=0 is limb)
ps plot (10^. }.taus);>(_100* }.p);(}.s) NB. plot log tau (exclude tau=0 point)