clear functions

net=readNet('dimerization.net');
net.species(1).initialAmount=2;
net.species(2).initialAmount=2;

X1values=0:7;
X2values=0:2;
%X1values=0:2;
%X2values=0:2;

[net1,expr1]=fsp(net,'X1',X1values,'X1_%d',[X1;X2]);
[net2,expr2]=fsp(net1,'X2',X2values,'X2_%d',expr1);

profile on
profile clear 

mdyn=momentDynamics(net2,2);

profile viewer

disp('Mu')
disp(mdyn.Mu.sym)
disp('Mu,dotMu')
disp([mdyn.Mu.sym,mdyn.dotMu.sym])
disp('A')
disp(mdyn.dotMu.A)
disp('B')
disp(mdyn.dotMu.B)
disp('barMu')
disp(mdyn.barMu.sym)

Tmax=1;
nMC=400;

mdyn.approxDotMu.sym=mdyn.dotMu.sym;
dynamics2mfile(net2,mdyn,'fun');
[t,mu]=ode23s(@(t,x)fun(x),[0,Tmax],mdyn.Mu.x0);

[handles,averages,stddevs]=plotCMoments(net2,mdyn,t,mu,expr2(1),'-',expr2(2),'--','X1_NaN','k-.','X2_NaN','k--',1);
legend(handles(1:3:end),'X1','X2','X1_{NaN}','X2_{NaN}');


