F={{0,1,0,0},{3 (w^2),0,0,2 w},{0,0,0,1},{0,-2 w,0,0}};
Q={{q1,0},{0,q2}};
G={{0,0},{1,0},{0,0},{0,1}};
P={{p11,p12,p13,p14},{p12,p22,p23,p24},{p13,p23,p33,p34},{p14,p24,p34,p44}};

(* Part (a): *)
H={{0,0,1,0}};

F.P + P.Transpose[F]+G . Q. Transpose[G]-(1/r)P.Transpose[H].H.P//MatrixForm

(* solving these is taking a very long time *)
Solve[{2 p12-((p11^2)/r)==0,p22-((p11 p12)/r)+2 p14 w+3 p11 (w^2)==0,p14+p23-((p11 p13)/r)==0,p24-((p11 p14)/r)-2 p12 w==0,
q1-((p23^2)/r)+4 p24 w+6 p12 (w^2)==0,p24-((p23 p33)/r)+2 p34 w+3 p13 (w^2)==0,-((p23 p34)/r)-2 p22 w+2 p44 w+3 p14 (w^2)==0,
2 p34-((p33^2)/r)==0,p44-((p33 p34)/r)-2 p23 w==0,
q2-((p34^2)/r)-4 p24 w==0
},{p11,p12,p13,p14,p22,p23,p24,p33,p34}]

(* create the block matrix in the solution of the Riccati equation *)
bigM=ArrayFlatten[{{-Transpose[F],(1/r)Transpose[H].H},{G.Q.Transpose[G],F}}];
MatrixForm[bigM]

bigPhi=Exp[bigM t];
MatrixForm[bigPhi]

bigPhi

(* extract out the pieces we need to solve for P(infty) *)
phiLambdaY={{1,1,1,1},{1,(E^(q1 t)),1,1},{1,1,1,1},{1,1,1,(E^(q2 t))}};
MatrixForm[phiLambdaY]
phiLambdaLambda={{1,(E^t),1,1},{(E^(3 t (w^2))),1,1,(E^(2 t w))},{1,1,1,(E^t)},{1,(E^(-2 t w)),1,1}};
MatrixForm[phiLambdaLambda]
phiLambdaY . phiLambdaLambda//Simplify//MatrixForm

(* this looks like it diverges as we take the limit of t going to infinity.  Can I use this method?  What are the eigenvalues of F? *)

Eigenvalues[F]

(* Since they do not have negative real parts we cannot us the above limit to determine P(infty) *)

(* Part (b): *)
H={{1,0,0,0}};

F.P + P.Transpose[F]+G . Q. Transpose[G]-(1/r)P.Transpose[H].H.P//MatrixForm