Possible to solve Algebraic Riccati Equation through ODE45? The Algebraic Riccati Equation is:
$$A'X + XA - XBR^{-1}B'X + Q = 0$$
And the Algebraic Riccati Difference Equation is:
$$A'X + XA - XBR^{-1}B'X + Q = \frac{dX}{dt}$$
Some times $\frac{dX}{dt}$ is just called $S$ in books. 
But is it not possible to solve the Algebraic Riccati Difference Equation through ODE45 and break the simulation when $-tol <= \frac{dX}{dt} <= tol$ ?
You may wonder "Why not use Schur's method to solve the Riccati equations?". Well, Schur's method is a very bad method and there is a reason why MATLAB and Octave are not using it. 
Sure, the method find the solution $X$ to the riccati equation, but build a control law $L$ for the LQR controller by using that solution is not what I recommending.
Right now I using fsolve, but that's a built in library in Octave and a function in Optimization package for MATLAB.
 A: Yes it is. Here is the solution:
Step 1: - Create the function. In this case it's CARE
function [value, isterminal, direction] = func(t, X, A, B, Q, R)

 X = reshape(X, size(A)); % Vector to matrix
 value = A'*X + X*A - X*B*inv(R)*B'*X + Q; % Values is the derivative of X
 value = value(:); % Turn value to a vector

 isterminal = 1;   % Stop the integration - An event!  
 direction  = 0;

Step 2: - Create your matrices
A =

     1     1
     2     1

>> B

B =

     1
     1

>> Q

Q =

     1     0
     0    1

>> R

R =

     5

Step 3: - Simulate
 >> [T X] = ode45(@(t,X)func(t, X, A, B, Q, R), [0 100], X0) ;

Step 4: - Find the last values for X and generate a matrix x
 >> x = reshape(X(size(X,1), :), size(A))

 x =

    8.7541    5.4230
    5.4230    5.0536

Step 5: - Check if the solution is almost zero
 >> A'*x + x*A - x*B*inv(R)*B'*x + Q

 ans =

     0.0022    0.0017
     0.0017    0.0013

Yes it is. 
Here is the plot how we can se that after a few seconds, the derivative of $X$ is almost zero.
I use this command
 plot(T, X)


....This is a knif....solution :)
A: 
You may wonder "Why not use Schur's method to solve the Riccati equations?". Well, Schur's method is a very bad method and there is a reason why MATLAB and Octave are not using it. 

Nope, that's precisely what matlab is using. It first checks whether $R$ is well-conditioned and if so forms the Hamiltonian and uses the QZ decomposition. If $R$ is not well-conditioned then it forms the Extended Hamiltonian and applies the Van Dooren's deflation. 
You can follow the logic from SciPy's implementation of solve_continuous_are which always uses the extended Hamiltonian instead of an $R$ dependent solution.
