You can solve the algebraic Riccati equation at each step by eigendecomposition of the associated Hamiltonian matrix (see for example https://en.wikipedia.org/wiki/Algebraic_Riccati_equation).

Within my Simulink model the Matlab Function block that implements the SDRE controller calls the function ARE_diag() with arguments the matrices composing the ARE. The function calls matrix_balance() to balance the Hamiltonian matrix before computing its spectral decomposition (the algorithm is taken from the book "Numerical Recipes": http://numerical.recipes/ ).

Otherwise, you can use iterative techniques like Kleinman algorithm (see http://ieeexplore.ieee.org/abstract/document/1024809/)

 @ Zeashan Khan 

Thank you very much for your info. In 'care' or 'dare' , the matrices A and B should be constant? Am I right?

If not, we cannot use care or dare to solve the problem?

Thank you. 

@ Giovanni Franzini

Thank you very much for your info. It's really helpful.

I didn't understand SDRE algorithm for nonlinear systems