function [PP,rms,tvec]=difflyap(A,B,C,M,P0,Q,nt,tmax);
% Solve differential Lyapunov equation,
% using double exponential
% 	inputs:
% A, B, C: state space operators
% M, P0: covariance of external excitations and initial condition
% Q: energy measure matrix
% nt, tmax: number of times and last time
% 	outputs:
% PP: state covariance matrix at times tvec
% rms: state rms at times tvec
% tvec: time vector

%%%% compute exponential matrix for dt
tvec=linspace(0,tmax,nt);
dt=tvec(2)-tvec(1); % smallest time interval
expmat=expm(A*dt); % state transition from 0 to dt

%%%% solve for steady state covariance
Pinf=mylyap(A,A',B*M*B');

%%%% march in time
n=size(A,1);
PP=zeros(n,n,nt); % for saving of covariance matrices
PP(:,:,1)=P0; % save the initial covariance
rms=zeros(size(C,1),nt); % for saving of the rms
rms(:,1)=sqrt(abs(diag(C*P0*C')));% rms of initial time
H=eye(n); % initial state transition matrix
for ind=2:nt;
  %%%% compute using double exponential
  H=H*expmat; % compute state transition matrix
  P=Pinf+H*(P0-Pinf)*H'; % march the covariance

  %%%% save results
  rms(:,ind)=sqrt(abs(diag(C*P*C')));
  PP(:,:,ind)=P;
  e=abs(trace(P*Q)); % mean kinetic energy at that time
  fprintf('%% it %3.0f / %3.0f %8.3f || %8.3f\n',ind,nt,tvec(ind),e);
end
