(*calculate the first k mometns of a PH distribution given with (α,T)*)
PHmoments[α_,T_,k_]:=
Module[{i,e,n,res},
n= Dimensions[T][[1]];
res={};
e=Table[1,{n},{1}];
Ti=Inverse[-T];
For[i=1,i≤k,i=i+1,
res=Append[res,i! (α.MatrixPower[Ti,i].e)[[1]] ];
];
res
]

(*  Calculates the canonical representation according to  
G.Horváth and M.Telek
On the canonical representation of phase type distributions 
Performance Evaluation,vol.66,no.8,pp.396-409, 2009
from any order 3 vector-matrix representation, if exists *)
PH3Canonical[v_,h0_] :=
Module[ {ii,λ,a0,a1,a2,a3,x1,x2,x3,x13,γu,γ0,γl,ϕ,f,p1,p2,p3,p,lower1,e},

If[Dimensions[h0][[1]] ≠3,Print["Dimension ≠ 3 !!! "];Return[] ];

e={{1},{1},{1}};
ii=IdentityMatrix[3];
λ=-Eigenvalues[h0];
Print["eigenvalues: ",λ];

If[(Re[λ[[3]]]<0 )||Abs[Im[λ[[3]]]]>0, 
 Print["eigenvalue error!!!"]; Return[{v,h0}]
];

a0=λ[[1]] λ[[2]] λ[[3]] ;
a1=λ[[1]] λ[[2]] + λ[[1]] λ[[3]] + λ[[2]] λ[[3]];
a2=λ[[1]] + λ[[2]] +  λ[[3]];
{a0,a1,a2}=Re[{a0,a1,a2}];

f[x_]:= x^3+a2  x^2+a1 x+a0;
γu=(a2+2 (Sqrt[a2^2-3 a1]))/3;
γ0=(a2+ (Sqrt[a2^2-3 a1]))/3;
If[λ[[1]]∈ Reals,γl=λ[[1]],γl=γ0];
lower1=- (((v.h0.h0.e)[[1]])/((v.h0.e)[[1]]));
Print["γu=",γu,", γl=",γl,", lower1=",lower1];
If[lower1>γl ,γl=lower1];

If[γl>γu ,Print["Max(γl,lower) > γu !!!!!!"];x1=γu,x1=γl];

x13=((-f[-x1])/(a0-f[-x1]))x1;
x2=(a2-x1+(Sqrt[(((a2-x1))^2)-4(x1^2 - a2 x1 +a1)]))/2;
x3=(a2-x1-(Sqrt[(((a2-x1))^2)-4(x1^2 - a2 x1 +a1)]))/2;
Print["x1=",x1,", x13=",x13,", x2=",x2,", x3=",x3];
{x13,x2,x3}=Re[{x13,x2,x3}];
(* Print["x1=",x1,", x13=",x13,", x2=",x2,", x3=",x3];*)
If[x13>x1,Print["x13 > x1 !!!!"] ];
If[x13<0,Print["x13 < 0 !!!!"] ];
p1=- h0.{{1},{1},{1}}/(x1-x13);p1=Transpose[p1][[1]];
p2=-(x1 ii + h0).h0.{{1},{1},{1}}/(x1-x13)/x2;p2=Transpose[p2][[1]];
p3=-(x2 ii + h0).(x1 ii + h0).h0.{{1},{1},{1}}/(x1-x13)/x2/x3;p3=Transpose[p3][[1]];
p=Transpose[{p1,p2,p3}];
Return[{v.p,{{-x1,0,x13},{x2,-x2,0},{0,x3,-x3}}}];
];

(* PH3 example with complex eigenvalues *) 
h0=(({{-3., 2, 0}, {0, -1, 1}, {1, 0, -2}}));
v={0.8, -0.3,0.5};
PHmoments[v,h0,5]
{α,A}=PH3Canonical[v,h0] ; {α,MatrixForm[A]}//Chop
PHmoments[α,A,5]

(* PH3 example with real eigenvalues and positive α_2 *) 
h0=(({{-3., 2, 1}, {0, -1, 1}, {1, 0, -2}}));
v={0.8, -0.4,0.4};
PHmoments[v,h0,5]
{α,A}=PH3Canonical[v,h0] ; {α,MatrixForm[A]}//Chop
PHmoments[α,A,5]

(* PH3 example with real eigenvalues and positive x_31 *) 
h0=(({{-3., 2, 1}, {0, -2, 2}, {1, 0, -2}}));
v={1.3, -1.0,0.7};
PHmoments[v,h0,5]
{α,A}=PH3Canonical[v,h0] ; {α,MatrixForm[A]}//Chop
PHmoments[α,A,5]

(* Non PH3 example 1 *) 
h0=(({{-3., 2, 1}, {0, -2, 2}, {1, 0, -2}}));
v={2.3, -2.0,0.7};
PHmoments[v,h0,5]
{α,A}=PH3Canonical[v,h0] ; {α,MatrixForm[A]}//Chop
PHmoments[α,A,5]

(* Non PH3 (non ME3) example 2 *) 
h0=(({{-3., 2, 1}, {0, -2, 2}, {-1, 0, -2}}));
v={1.3, -1.0,0.7};
PHmoments[v,h0,5]
{α,A}=PH3Canonical[v,h0] ; {α,MatrixForm[A]}//Chop
PHmoments[α,A,5]
Plot[{α.MatrixExp[A x].(-A.Table[{1},{Length[α]}])},{x,0,10},AxesOrigin->{0,0}]

(* Valid ME3, but non PH3 example 3 *) 
h0=(({{-1, 1, 0}, {0, -2, 2}, {2.5, 0, -3}}));v={0.8,-0.4,0.6};
PHmoments[v,h0,5]
{α,A}=PH3Canonical[v,h0] ; {α,MatrixForm[A]}//Chop
PHmoments[α,A,5]
Plot[{α.MatrixExp[A x].(-A.Table[{1},{Length[α]}])},{x,0,10},AxesOrigin->{0,0}]