nonsymmetric differential Riccati matrix equations

The resolution of the nonsymmetric differential matrix Riccati equations
%
% dY(t)/dt = AY + YB - YCY + Q (*)
% Y(t0) = Y0
%
% using backward differentiation formula method.
%
% Input: A : size matrix (n,n).
% B : size matrix (p,p).
% C : size matrix (p,n).
% D : size matrix (n,p).
% Y0 : size matrix (n,p).
% tf , t0 : temps.
% Output:
% Ytf : solution to the inslant tf of (*)
% Y : solution to (*) for t\in[t0 tf]
% t
% Author : LAKHLIFA SADEK.
% Last modification : 29/09/2019
% E-mail: lakhlifasdek@gmail.com; sadek.l@ucd.ac.ma
% ORCID : https://orcid.org/0000-0001-9780-2592
%
%
% Example
%
% t0=0, tf=1
%
% A =
% 1 2
% 3 5
%
% B =
% 2 2 1
% 2 3 2
% 2 1 1
%
% C =
% 1 4
% 1 3
% 2 3
%
% D =
% 1 5 3
% 2 6 2
%
%
% Y0 =
%
% 1 4 3
% 2 3 -1
%
% [t,Y,Ytf] = BDF_Diff_Riccati(A,B,C,D,Y0,t0,tf)
%
% Y(tf) solution
%
% Ytf =
%
% 0.0790 0.3842 1.3520
% 1.2402 1.8543 0.7197
%
% for tracer the solution Y_{1,1}(t) to fonctuion t
% plot(t,Y(1,:))
% xlabel('Temps t\in[0 1]')
% ylabel('The solution Y_{1,1}(t)')
% legend('BDF method').