% This example is discussed in the documentation of StochDynTools

clear functions
clear all

%%%%%%%%%%%%%%%%
% read .net file
%%%%%%%%%%%%%%%%
net=readNet('enzyme_with_Et_St.net');

%%%%%%%%%%%%%%%%%%%%%%
% 2nd order trunaction
%%%%%%%%%%%%%%%%%%%%%%
[mdyn2,x20,texrules]=closureDynamics(net,2,'derMatch','fun');

%%%%%%%%%%%%
% simulation 
%%%%%%%%%%%%
Tmax=300;                      % maximum simulation time
[t,mu]=ode23s(@(t,x)fun(x),[0,Tmax],x20);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot of uncentered moments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1)
plot(t,mu)
legend(fun(),'location','best','interpreter','none')
print -depsc uncentered.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot of centered moments
%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2)
h=plotCMoments(net,mdyn2,t,mu,'E','b-','St','g-','Et','k-');
legend(h(1:3:end),'E[E]\pm{}Std[E]','E[St]\pm{}Std[St]','E[Et]','location','best');
print -depsc centered.eps

%%%%%%%%%%%%%%%%%%%%%%
% plot of distribution
%%%%%%%%%%%%%%%%%%%%%%
figure(3)
x=0:.01:4;
pdf=getDistribution(net,mdyn2,mu(end,:),'St',x,'lognormal');
plot(x,pdf)
print -depsc pdf.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1st order trunaction & latex output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
[mdyn1,x10,texrules]=closureDynamics(net,1,'derMatch');
dotMu=expand(mdyn1.dotMu.approxDotMu);
% dotMu=expand(subsParameters(net,dotMu));
latexmacros=mylatex()
latexMu=mylatex(mdyn1.Mu.sym,texrules)
latexDotMu=mylatex(dotMu,texrules)

%%%%%%%%%%%%%%%%%%%
% equilibrium point
%%%%%%%%%%%%%%%%%%%
equations=cell(mdyn1.dotMu.approxDotMu);
variables=cell(mdyn1.Mu.sym);
solution=solve(equations{:},variables{:})



