Procedure to match mean, squared coefficient of variation, the third order moment expression l3 (optional), and the lag-1 autocorrelation.

MAPMatch[m1_,c2_,l3_,ρ1_,rr_:None]:=Module[{r,m11,m12,p1,p2,cv21,cv22,l31,l32,α1,A1,α2,A2,N1,N2,h1,h2,D0,D1,m21,m22,m31,m32},
If[ρ1≥0,
If[rr=!=None && rr≤0 && rr≥1,Print["Parameter r is out of range!"];Abort[]];
If[rr===None, 
r=((2 ρ1)/(1+ρ1));
p1= (1/(1+c2))(1-((1+ρ1)/2));
p2= (c2/(1+c2))(1-((1+ρ1)/2));
,
r=rr;
p1= (1/(1+c2))(1-((ρ1)/r));
p2= (c2/(1+c2))(1-((ρ1)/r));
];
m11 = m1(1-(Sqrt[r]));
m12 = m1(1+c2 (Sqrt[r]));
If[l3===None,
cv21=(((Sqrt[c2]) (1+c2) (1+(Sqrt[r])))/(1-(Sqrt[c2]) (-1+(Sqrt[r]))+c2 (Sqrt[r])));
cv22=-((c2 (1+c2) (-1+r))/((1+c2 (Sqrt[r])) (1-(Sqrt[c2]) (-1+(Sqrt[r]))+c2 (Sqrt[r]))));
(*Print[r,", ",m11,", ", m12, ", ", cv21, ", ",cv22];*)
{α1,A1}=PHFrom2Moments[m11,cv21];
{α2,A2}=PHFrom2Moments[m12,cv22];
,
cv21=((c2+(Sqrt[r]))/(1-(Sqrt[r])));
cv22=c2((1-(Sqrt[r]))/(1+c2(Sqrt[r])));
l31=(((1+c2) l3)/(c2(1-(Sqrt[r]))+(Sqrt[(1 +c2 (Sqrt[r]))c2(1-(Sqrt[r]))])));
l32=(((1+c2) l3)/((1+c2(Sqrt[r]))+(Sqrt[(1 +c2 (Sqrt[r]))c2(1-(Sqrt[r]))])));
(*Print[r,", ",m11,", ", m12, ", ", cv21, ", ",cv22,", ",l31,", ", l32];*)
If[cv21≤0||cv22≤0||l31≤0 || l32≤0, Print["Not feasible. c2=",c2,", l3=",l3, ", r=",r]; Return[{{},{}}]];
m21 = (cv21+1)(m11^2);
m22 = (cv22+1)(m12^2);
m31 = (l31+1)(m21^2)/m11;
m32 = (l32+1)(m22^2)/m12;
{α1,A1}=PHFrom3Moments[m11,m21,m31];
{α2,A2}=PHFrom3Moments[m12,m22,m32];
],
If[c2≥1,
If[rr=!=None && rr≤0 && rr≥1/c2,Print["Parameter r is out of range!"];Abort[]];
If[rr===None, 
r=-((2 ρ1)/(1-c2 ρ1));
p1= (1/2)(1+((1-c2 ρ1)/2));
p2= (1/2)(1+((1-c2 ρ1)/2));
,
r=rr;
p1= (1/2)(1-((ρ1)/r));
p2= (1/2)(1-((ρ1)/r));
];
,
If[rr=!=None && rr≤0 && rr≥1,Print["Parameter r is out of range!"];Abort[]];
If[rr===None, 
r=-((2 ρ1)/(1-ρ1));
p1= (1/2)(1+((1-ρ1)/2));
p2= (1/2)(1+((1-ρ1)/2));
,
r=rr;
p1= (1/2)(1-((ρ1)/r));
p2= (1/2)(1-((ρ1)/r));];
];
m11 = m1(1-(Sqrt[c2 r]));
m12 = m1(1+ (Sqrt[c2 r]));
If[l3===None,
cv21=c2 ((1-r)/(1-(Sqrt[c2 r])));
cv22=c2 ((1-r)/(1+(Sqrt[c2 r])));
(*Print[r,", ",m11,", ", m12, ", ", cv21, ", ",cv22];*)
{α1,A1}=PHFrom2Moments[m11,cv21];
{α2,A2}=PHFrom2Moments[m12,cv22];
,
cv21=((c2+(Sqrt[c2 r]))/(1-(Sqrt[c2 r])));
cv22=((c2-(Sqrt[c2 r]))/(1+(Sqrt[c2 r])));
l31=2 l3 (1/(1-(Sqrt[c2 r])+(Sqrt[1-c2 r])));
l32=2 l3 (1/(1+(Sqrt[c2 r])+(Sqrt[1-c2 r])));
If[cv21≤0||cv22≤0||l31≤0 || l32≤0, Print["Not feasible. c2=",c2,", l3=",l3, ", r=",r]; Return[{{},{}}]];
(*Print[r,", ",m11,", ", m12, ", ", cv21, ", ",cv22,", ",NumberForm[l31,50],", ", l32];*)
m21 = (cv21+1)(m11^2);
m22 = (cv22+1)(m12^2);
m31 = (l31+1)(m21^2)/m11;
m32 = (l32+1)(m22^2)/m12;
{α1,A1}=PHFrom3Moments[m11,m21,m31];
{α2,A2}=PHFrom3Moments[m12,m22,m32];
];
];
N1 = Length[α1];
N2 = Length[α2];
h1=Table[{1},{N1}];
h2=Table[{1},{N2}];
D0=Table[0,{N1+N2},{N1+N2}];
D0[[1;;N1,1;;N1]]=A1;
D0[[N1+1;;N1+N2,N1+1;;N1+N2]]=A2;
D1=Table[0,{N1+N2},{N1+N2}];
D1[[1;;N1,1;;N1]]=-A1.h1.{α1}(1-p1);
D1[[1;;N1,N1+1;;N1+N2]]=-A1.h1.{α2} p1;
D1[[N1+1;;N1+N2,1;;N1]]=-A2.h2.{α1} p2;
D1[[N1+1;;N1+N2,N1+1;;N1+N2]]=-A2.h2.{α2}(1-p2);
Return[{D0,D1}];
];

Procedure to construct a PH distribution based on the mean and the squared coefficient of variation

PHFrom2Moments[m_,cv2_]:=Module[{N,A,p,α,λ,i},
λ=1/m;
N=Max[Ceiling[1/cv2],2];
p=(1/(cv2+1+((cv2-1)/(N-1))));
A=DiagonalMatrix[Table[-λ p N,{N}]];
For[i=1,i<N,i++,A[[i,i+1]]=-A[[i,i]]];
A[[N,N]]=-λ N;
α=Table[0,{N}];
α[[1]]=p;
α[[N]]=1-p;
Return[{α,A}];
];

Procedure to construct a PH distribution based on 3 moments

PHFrom3Moments[m1_,m2_,m3_]:=Module[{mym1,n2,n3,αm,A},
	If[Det[({ {1, m1}, {m1, m2}})]<=0,
	Print["Non-positive variance"]; 
	Return[{{0},{{0}}}]];
	If[Det[({ {m1, m2}, {m2, m3}})]<=0,
	Print["No positive distribution with such moments"]; 
	Return[{{0},{{0}}}]];
	{mym1,n2,n3}=MomsToNormMoms[{m1,m2,m3}];
	{{αm},A}=MOM3Construct[mym1,n2,n3];
	Return[{αm,A}];
]

MomsToNormMoms[ moms_] :=Module[ {isize, i, normmoms ,m},
	isize=Dimensions[moms][[1]];
	normmoms={moms[[1]]};
	For[ i= 2,i<= isize,
	m=moms[[i]]/(moms[[i-1]] moms[[1]]);
	normmoms=Append[normmoms,m];i++];
	Return[normmoms];
];

BuToolsReal[x_,n_]:=If[Im[x]==0||Abs[Re[x]/Im[x]]>10^n,Return[Re[x]],Return[x]];

(* Construct minimal PH dist with normalized moments n1,n2,n3  *)
MOM3Construct[n1_,n2_,n3_]:=Module[{n,b,b1,b2,a1,a2,c0,c1,c2,c3,c4,K1,K2,K3,K4,K5,K6,P1,P2,P3,P41,P42,m1,m2,m3,p1,p2,lambda1,lambda2,mu1,mu2,result,fs,K7,K8,K9,i1,ff,a,p,lambda,mu},
n=2;
While[(n+1)/n>n2 || MOM3lowerbound[n,n2]≥n3||MOM3upperbound[n,n2]<n3,n=n+1];
If[n3<=2n2-1,
(*Print["Below the Erlang-line:    ErlangExp"];*)
	b=((2(4-n(3 n2-4)))/(n2(4+n-n n3)+(Sqrt[n n2])(Sqrt[12 (n2^2)(n+1)+16 n3(n+1)+n2(n(n3-15)(n3+1)-8(n3+3))])));
a=(( b n2-2)(n-1)b)/((b-1)n);
	p=(b-1)/a;
	lambda=(p a+1)/n1;
	mu=p/(n1-1/lambda) (n-1);
	If[0<= p<= 1&&lambda>0&&mu>0,Return[MOM3ErlangExp[n,p,lambda,mu]]];
	result={n,p,lambda,mu};
	,
	(*Print["Above the Erlang-line:    ExpErlang"];*)
	c4=n2 (3n2-2n3)(n-1)^2;
	c3=2n2 (n3-3)(n-1)^2;
	c2=6(n-1)(n-n2);
	c1=4n (2-n);
	c0=n (n-2);
fs=NSolve[c4 f^4+c3 f^3+c2 f^2+c1 f+c0==0,f];
result={n};
For[i1=1,i1<=4,i1=i1+1,
ff=fs[[i1,1,2]];
If[(n-1)(n2 ff^2-2ff+2)-n==0,Continue[]];
a=(2(ff-1)(n-1))/((n-1)(n2 ff^2-2ff+2)-n);
p=N[(ff-1)a,Precision[(ff-1)a]];
lambda=N[(a+p)/n1,Precision[(a+p)/n1]];
mu=N[(n-1)/(n1-p/lambda),Precision[(n-1)/(n1-p/lambda)]];
result=Append[result,p];
result=Append[result,lambda];
result=Append[result,mu];
p=BuToolsReal[p,Precision[p]3/4];
lambda=BuToolsReal[lambda,Precision[lambda]3/4];
mu=BuToolsReal[mu,Precision[mu]3/4];
If[Element[p,Reals]&&Element[lambda,Reals]&&Element[mu,Reals],
If[0<= p<= 1&&lambda>0&&mu>0,Return[MOM3ExpErlang[n,p,lambda,mu]]];
];
];
];
Print["Unexpected termination of MOM3Construct!"];
Return[result];
];

MOM3ExpErlang[n_,p1_,lambda1_,mu1_]:= Module[{AA,TT,resTT,FF,m1,m2,m3,i},
TT=-mu1 IdentityMatrix[n];
For[i=2,i<n,i=i+1,TT[[i,i+1]]=mu1];
TT[[1,1]]=-lambda1;
TT[[1,2]]=lambda1;
AA=Table[0,{i,1,1},{j,1,n}];
AA[[1,1]]=p1;
AA[[1,2]]=1-p1;
Return[{AA,TT}];
];

MOM3ErlangExp[n_,p1_,lambda1_,mu1_]:= Module[{AA,TT,resTT,FF,m1,m2,m3},
TT=-mu1 IdentityMatrix[n];
For[i=1,i<n,i=i+1,TT[[i,i+1]]=mu1];
TT[[n,n]]=-lambda1;
AA=Table[0,{i,1,1},{j,1,n}];
AA[[1,1]]=p1;
AA[[1,n]]=1-p1;
Return[{AA,TT}];
];

MOM3upperbound[n_,n2_]:=Module[{ub,u},
ub=Which[
n2<(n+1)/n,-Infinity,
n2<= n/(n-1),
	u=1/(n^2 n2) (2(n-2)(n n2-n-1)Sqrt[1+(n (n2-2))/(n-1)]+(n+2)(3n n2-2n-2));
	u,
True,Infinity]];

MOM3lowerbound[n_,n2_]:=Module[{lb,l,p,a},
lb=Which[
n2<(n+1)/n,Infinity,
n2< (n+4)/(n+1),
	p=((n+1)(n2-2))/(3n2 (n-1)) ((-2Sqrt[n+1])/Sqrt[-3n n2+4n+4]-1);
	a=(n2-2)/(p (1-n2)+Sqrt[p^2+(p n (n2-2))/(n-1)]);
	l=((3+a)(n-1)+2a)/((n-1)(1+a p))-(2a (n+1))/(2(n-1)+a p (n a+2n-2));
	l,
True,(n+1)/n n2]];

Test of the procedure:

(*{D0,D1}=MAPMatch[1.1`50, 4`50, 3`50,-0.17`50];*)
{D0,D1}=MAPMatch[1.2`50, 0.6`50, 8.1`50,-0.2`50];
(*{D0,D1}=MAPMatch[1.2`50, 0.4`50, 0.01`50,-0.20`50];*)
Dimensions[D0]
P=Inverse[-D0].D1;
α=Table[(Subscript[x, i]),{i,1,Dimensions[D0][[1]]}];
h=Table[{1},{i,1,Dimensions[D0][[1]]}];
α=α/.Solve[{α.P==α,α.h==1},α];
(* Check m1 value: *)
α.Inverse[-D0].h
(* Check c2 value: *)
2 α.MatrixPower[Inverse[-D0],2].h/(((α.Inverse[-D0].h))^2)-1
(* Check l3 value: *)
(6 α.MatrixPower[Inverse[-D0],3].h)(α.Inverse[-D0].h)/(((2 α.MatrixPower[Inverse[-D0],2].h))^2)-1
(* Check ρ1 value: *)
((α.Inverse[-D0].P.Inverse[-D0].h - (((α.Inverse[-D0].h))^2))/(2 α.MatrixPower[Inverse[-D0],2].h - (((α.Inverse[-D0].h))^2)))