clear functions
Tmax=5;
net = readNet('mather.net');
mdyn = closureDynamics(net,2,'dm');

mdyn.Mu.sym
mdyn.dotMu.sym

muEq = momentEquilibrium(net,mdyn)

X1minmax=[0,20];
X2minmax=[0,20];
nMC=10000;

figure(1);clf
[exact,approx,MC]=plotCompare(net,2,5,Tmax,...
    {{'derivative matching',X1minmax},{'MCsample',1,'X1','c-'},{'dm',2,'X1','b-'}},...
    {{'zero cumulants',X1minmax},{'zc',1,'X1','k-'},{'zc',2,'X1','b-'}},...
    {{'low dispersion',X1minmax},{'ld',1,'X1','k-'},{'ld',2,'X1','b-'}},...
    {{'quasi deterministic',X1minmax},{'qd',1,'X1','k-'},{'qd',2,'X1','b-'}},...
    {{'distribution (Monte Carlo)',X1minmax},'MCdistribution',nMC,'X1'},...
    {{'derivative matching',X2minmax},{'MCsample',1,'X2','c-'},{'dm',2,'X2','b-'}},...
    {{'zero cumulants',X2minmax},{'zc',1,'X2','k-'},{'zc',2,'X2','b-'}},...
    {{'low dispersion',X2minmax},{'ld',1,'X2','k-'},{'ld',2,'X2','b-'}},...
    {{'quasi deterministic',X2minmax},{'qd',1,'X2','k-'},{'qd',2,'X2','b-'}},...
    {{'distribution (Monte Carlo)',X2minmax},'MCdistribution',nMC,'X2'});


