(* solve for the steady state uncertainty pInfinity *)
pInfRoots=Solve[((-(b^2))/r)(p^2)+2 a p+q==0,p]

(* we must take the postive root since p_{Infty} must be postive *)
pInf=pInfRoots[[2,1,2]]

(* the transformation z = 1/(p-pInf) gives a coefficient of z[t] given by *)
zCoef=-Simplify[2 a - 2 ((b^2)/r)pInf]

(* solve for z[t]  *)
zSol=DSolve[{z'[t]==2 beta z[t]+((b^2)/r),z[0]==z0},z[t],t][[1,1,2]]

(* from z[t] compute what the solution for p[t] is *) 
Clear[pInf];
pSol=Simplify[pInf+(1/zSol)]

(* find out what the value of z0 is in terms of p0 the initial uncertainty in the state *)
p0Eval=pSol/.{t->0}
z0InTermsOfP0=Solve[p0Eval==p0,z0]//Simplify

(* put this value back into the solution for p[t] *)
pSol2=pSol/.z0InTermsOfP0//Simplify

(* convert to hyperbolic functions *)
pSol3=ExpToTrig[pSol2[[1]]]//Together//Simplify

(* put in pInf in terms of beta and simlify *)
pSol3/.{pInf->(r/(b^2))(beta + a)}//Simplify

(* this last expression can be manipulated into the form claimed in the book *)