kFullStep[pPlus0_]:=Module[{pMinus,pPlus1Inv,pPlus1},
(* this does the algebra for one step of the discrete Kalman filter starting from p(+), producing the next p(-), and the next p(+) *)
(* pMinus=Simplify[pPlus0+q Delta];*)
 pMinus=Simplify[pPlus0+g r0];(* gamma r0 *)
pPlus1Inv=Together[(1/pMinus)+(1/r0)];
pPlus1=(1/pPlus1Inv);
Return[{pMinus,pPlus1}]; (* returns p(-) and p(+) *)
];

(* Lets do the two steps of the forward solution *)
p0Plus=(1/smallEps);
p1=kFullStep[p0Plus]
p2=kFullStep[p1[[2]]]

(* what do these expressions evaluate when epsilon is zero? *)
p1/.{smallEps->0}
p2/.{smallEps->0}

We will do the backwards smoothing by hand now when epsilon is not exactly equal to zero.

(* compute A1 first *)
A1=Simplify[((p1[[2]])/(p2[[1]]))] (* this is p1(+)/p2(-) *)
A1/.{smallEps->0}

(* now compute P_{1|2} *) 
p1Given2=Simplify[p1[[2]]+A1 ( p2[[2]]-p2[[1]] ) A1]
p1Given2/.{smallEps->0}

(* compute A0 next *)
A0=Simplify[(p0Plus/(p1[[1]]))] (* this is p0(+)/p1(-) *)
A0/.{smallEps->0}

(* now compute P_{0|2} *)
p0Given2=Simplify[p0Plus+A0 (p1Given2 - p1[[1]]) A0]
p0Given2/.{smallEps->0}