Your Little CMB code by Enrico Pajer
(based on Seljak’s [astro-ph/9406050], on Zaldarriaga’s and Baumann’s exercise, and corrections suggested by Essex)

Definitions

ClearAll[lsqCl,ΔRsq,y,η,τfct]
(*physical fractions ω = Ω (h^2) at z=0*)
(Subscript[ω, m])= 0.13; (Subscript[ω, b]) = 0.02; (Subscript[ω, c])=(Subscript[ω, m])-(Subscript[ω, b]); ωγ=2.47 (10^(-5)); ωΛ=0.728 (((0.7))^2); 
(Subscript[a, dec])=(1100^(-1));  (*redshift at photon decoupling*)
(Subscript[a, eq])= (((4.2*(10^(-5))))/(Subscript[ω, m])); (*redshift at matter-radiation equality, with neutrinos!*)
(Subscript[f, ν]) = 1.68;     (*total radiation density  (Subscript[ρ, r])=(Subscript[ρ, γ])+(Subscript[ρ, ν])=[1+3(7/8)((((4/11)))^(4/3))](Subscript[ρ, γ])=(Subscript[f, ν]) (Subscript[ρ, γ])*)

τfct[af_]:=(1/Mpc)(((100 (km/sec)(1/Mpc)))^(-1))(300 (10^3) (km/sec)) NIntegrate[(1/((a^2) (Sqrt[(a^(-4))ωγ (Subscript[f, ν])+(a^(-3))((Subscript[ω, m]))+ωΛ]))),{a,0,af}];
(Subscript[τ, L])=τfct[1]-τfct[(Subscript[a, dec])](*distance to the surface of last scattering in Mpc*)

(*Rescaled variables*)
α= (Sqrt[((Subscript[a, dec])/(Subscript[a, eq]))]);
(*conformal time at recombination in a would-be matter dominated universe in units of Mpc*)
(Subscript[τ, r])= (2.997 (10^3))*(Sqrt[4(Subscript[a, dec])/((Subscript[ω, m]))])   ;
(Subscript[x, i])= (10^(-5));  (*starting time of the code in dimensionless units x=(τ/(Subscript[τ, r])), corresponding to z=(10^8)*)
(Subscript[x, dec])= ((Sqrt[(α^2)+1])-1)/α ;(*dimensionless time at decoupling*)
(Subscript[κ, min])=0.01;   (*dimensionless wavenumber k=κ/(Subscript[τ, r])=(10^(-5)), largest scales*)
(Subscript[κ, max ])=300;   (*shortest scales computed by the code, k=κ/(Subscript[τ, r])=0.6. The largest argument I will use is κ=(x l (Subscript[τ, r])/(Subscript[τ, L]))<257*)

(*Damping*)
σ= 0.03; (*width of visibility function, see astro-ph/9406050*)
(Subscript[x, s])=0.6((Subscript[ω, m])^0.25)((Subscript[ω, b])^(-0.5))((Subscript[a, dec])^0.75);(*Silk damping scale*)
(Subscript[k, d])= ((((2((Subscript[x, s])^2) + (σ^2)((Subscript[x, dec])^2))((Subscript[τ, r])^2)))^(-0.5)) ;(*Final damping scale*)

Integration

y[x_]:=2α x+(α^2)(x^2);  (*analytical background solution*)
η[x_]:=2α(1+α x)/(2α x+(α^2)(x^2)); (*rescaled Hubble analytical soltion*)
(*Rescaled equations and initial conditions*)
ICs= {ϕ[(Subscript[x, i])]==1,(Subscript[δ, γ])[(Subscript[x, i])]==-2 ϕ[(Subscript[x, i])],(Subscript[δ, c])[(Subscript[x, i])]==(3/4)(Subscript[δ, γ])[(Subscript[x, i])],(Subscript[v, c])[(Subscript[x, i])]==(Subscript[v, γ])[(Subscript[x, i])],(Subscript[v, γ])[(Subscript[x, i])]==-(1/4)((κ )/(η[(Subscript[x, i])]))(Subscript[δ, γ])[(Subscript[x, i])]};
equations= {(Subscript[δ, c])'[x]==-κ (Subscript[v, c])[x] + 3ϕ'[x],
(Subscript[v, c])'[x]==-η[x] (Subscript[v, c])[x] + κ ϕ[x],
(Subscript[δ, γ])'[x]==-(4/3)κ (Subscript[v, γ])[x] + 4 ϕ'[x],
(Subscript[v, γ])'[x]== κ ϕ[x]+(-0.75((Subscript[ω, b])/(Subscript[ω, m])) (Subscript[f, ν]) y[x]η[x](Subscript[v, γ])[x] + (1/4)κ (Subscript[δ, γ])[x])(((1+0.75((Subscript[ω, b])/(Subscript[ω, m])) (Subscript[f, ν]) y[x]))^(-1)),
ϕ'[x]==-η[x]ϕ[x]+ 1.5(((η[x])^2)/κ)((4/3 + ((Subscript[ω, b])/(Subscript[ω, m])) y[x])(Subscript[v, γ])[x] + ((Subscript[ω, c])/(Subscript[ω, m]))y[x](Subscript[v, c])[x])/(1+y[x]),
ICs};
(*Integration module*)
integrate:=
Flatten[Module[{},
sol=NDSolve[equations,{(Subscript[δ, c]),(Subscript[δ, γ]),(Subscript[v, c]),(Subscript[v, γ]),ϕ},{x,(Subscript[x, i]),(Subscript[x, dec]) }];
     {ϕ[(Subscript[x, dec])]+(Subscript[δ, γ])[(Subscript[x, dec])]/4/.sol ,(Subscript[v, γ])[(Subscript[x, dec])]/.sol}]]

(*Compute the monopole and dipole at decoupling*)
(Subscript[Θ, 0])=Interpolation[Table[{κ,integrate[[1]]},{κ,(Subscript[κ, min]),(Subscript[κ, max]),1}]];
(Subscript[Θ, 1])=Interpolation[Table[{κ,-integrate[[2]]},{κ,(Subscript[κ, min]),(Subscript[κ, max]),1}]];

LogLinearPlot[{(Subscript[Θ, 0])[k (Subscript[τ, r])],(Subscript[Θ, 1])[k (Subscript[τ, r])]/(Sqrt[3])},{k,(Subscript[κ, min])/(Subscript[τ, r]),(Subscript[κ, max])/((Subscript[τ, r]) 2)},PlotStyle->{Red,Blue},GridLines->{{0},{0}},AxesLabel->{"k[(Mpc^(-1))]","(Subscript[Θ, 0]),(Subscript[Θ, 1])"},PlotLabel->"Monopole and dipole at decoupling"]

(*Approximations of the spherical Bessel functions (Subscript[j, l])*)
w1[x_]:= (1/(2 (x(Sqrt[(x^2)-1]))))   (*(l^2) ((Subscript[j, l])^2)[x l]*)
w2[x_]:= w1[x]((x^2)-1)/(x^2)   (*(l^2) ((((Subscript[j, l])'[x l]))^2)*)
;

(*Primordial power spectrum*)
ΔSqRpivot=2.2 (10^(-9));ns=0.967;kpivot=0.05;
ΔRsq[k_]:=ΔSqRpivot ((((k/kpivot)))^(ns-1))

(*Projecting the power spectrum at decoupling onto the Cl's*)
lsqCl[l_]:= NIntegrate[((((Subscript[Θ, 0])[x l (Subscript[τ, r])/(Subscript[τ, L])])^2)w1[x]+ (((Subscript[Θ, 1])[x l (Subscript[τ, r])/(Subscript[τ, L])])^2)w2[x])*Exp[-(((x l (Subscript[τ, r])/(Subscript[τ, L])))^2)(2((Subscript[x, s])^2) + (σ^2)((Subscript[x, dec])^2))]((ΔRsq[x l/(Subscript[τ, L])])/x),{x,1,5}];
cmb=Interpolation[Table[{l,lsqCl[l]},{l,1,1500,20}]];

(*...and finally, your own Cl's!*)
Plot[cmb[l],{l,1,1500},AxesLabel->{"l","l(l+1)(Subscript[C, l])"}]