2 added 3 characters in body edited Oct 20 '17 at 22:04 OkkesDulgerci 6,06711 gold badge1111 silver badges2121 bronze badges Here is improved version  exact[\[Lambda]_, \[Mu]_, initialPopulation_] := Module[{}, DSolveValue[{p'[t] == (\[Lambda] - \[Mu]) p[t], p[0] == initialPopulation}, p[t], t]]; birthDeath[\[Lambda]_, \[Mu]_, initialPopulation_, numOfReaction_] := NestList[(\[CapitalDelta]t1 = RandomVariate[ExponentialDistribution[\[Lambda] #[[2]]]]; \[CapitalDelta]t2 = RandomVariate[ExponentialDistribution[\[Mu] #[[2]]]]; {\[CapitalDelta]t = Min[\[CapitalDelta]t1, \[CapitalDelta]t2]\[CapitalDelta]t2]; {#[[1]] + \[CapitalDelta]t, If[\[CapitalDelta]t1 < \[CapitalDelta]t2, #[[2]] + 1, #[[2]] - 1]}) &, {0, initialPopulation}, numOfReaction]   With[{\[Lambda] = 3, \[Mu] = 1, initialPopulation = 10, numOfReaction = 2001000, numOfsim = 10}, sim = Table[ birthDeath[\[Lambda], \[Mu], initialPopulation, numOfReaction], numOfsim]; sim[[All, All, 1]] = Accumulate /@ sim[[All, All, 1]]; Show[ListStepPlot[sim, PlotLegends -> {"Simulation"}, PlotStyle -> Directive[AbsoluteThickness[0.2]], Frame -> True, PlotTheme -> "Detailed", FrameLabel -> {"Time", "Population"}, ImageSize -> Large], Plot[Evaluate@exact[\[Lambda], \[Mu], initialPopulation], {t, 0, Max@sim[[All, -1, 1]]}, PlotStyle -> {Black, Thick}, PlotLegends -> {"ODE"}], ListLinePlot[Mean@sim, PlotLegends -> {"Avg of Simulation"}, PlotStyle -> {Red, Dashed}]]]  Here is improved version exact[\[Lambda]_, \[Mu]_, initialPopulation_] := Module[{}, DSolveValue[{p'[t] == (\[Lambda] - \[Mu]) p[t], p[0] == initialPopulation}, p[t], t]]; birthDeath[\[Lambda]_, \[Mu]_, initialPopulation_, numOfReaction_] := NestList[(\[CapitalDelta]t1 = RandomVariate[ExponentialDistribution[\[Lambda] #[[2]]]]; \[CapitalDelta]t2 = RandomVariate[ExponentialDistribution[\[Mu] #[[2]]]]; {Min[\[CapitalDelta]t1, \[CapitalDelta]t2], If[\[CapitalDelta]t1 < \[CapitalDelta]t2, #[[2]] + 1, #[[2]] - 1]}) &, {0, initialPopulation}, numOfReaction]   With[{\[Lambda] = 3, \[Mu] = 1, initialPopulation = 10, numOfReaction = 200, numOfsim = 10}, sim = Table[ birthDeath[\[Lambda], \[Mu], initialPopulation, numOfReaction], numOfsim]; sim[[All, All, 1]] = Accumulate /@ sim[[All, All, 1]]; Show[ListStepPlot[sim, PlotLegends -> {"Simulation"}, PlotStyle -> Directive[AbsoluteThickness[0.2]], Frame -> True, PlotTheme -> "Detailed", FrameLabel -> {"Time", "Population"}, ImageSize -> Large], Plot[Evaluate@exact[\[Lambda], \[Mu], initialPopulation], {t, 0, Max@sim[[All, -1, 1]]}, PlotStyle -> {Black, Thick}, PlotLegends -> {"ODE"}], ListLinePlot[Mean@sim, PlotLegends -> {"Avg of Simulation"}, PlotStyle -> {Red, Dashed}]]]  Here is improved version  exact[\[Lambda]_, \[Mu]_, initialPopulation_] := Module[{}, DSolveValue[{p'[t] == (\[Lambda] - \[Mu]) p[t], p[0] == initialPopulation}, p[t], t]]; birthDeath[\[Lambda]_, \[Mu]_, initialPopulation_, numOfReaction_] := NestList[(\[CapitalDelta]t1 = RandomVariate[ExponentialDistribution[\[Lambda] #[[2]]]]; \[CapitalDelta]t2 = RandomVariate[ExponentialDistribution[\[Mu] #[[2]]]]; \[CapitalDelta]t = Min[\[CapitalDelta]t1, \[CapitalDelta]t2]; {#[[1]] + \[CapitalDelta]t, If[\[CapitalDelta]t1 < \[CapitalDelta]t2, #[[2]] + 1, #[[2]] - 1]}) &, {0, initialPopulation}, numOfReaction] With[{\[Lambda] = 3, \[Mu] = 1, initialPopulation = 10, numOfReaction = 1000, numOfsim = 10}, sim = Table[ birthDeath[\[Lambda], \[Mu], initialPopulation, numOfReaction], numOfsim]; Show[ListStepPlot[sim, PlotLegends -> {"Simulation"}, PlotStyle -> Directive[AbsoluteThickness[0.2]], Frame -> True, PlotTheme -> "Detailed", FrameLabel -> {"Time", "Population"}, ImageSize -> Large], Plot[Evaluate@exact[\[Lambda], \[Mu], initialPopulation], {t, 0, Max@sim[[All, -1, 1]]}, PlotStyle -> {Black, Thick}, PlotLegends -> {"ODE"}], ListLinePlot[Mean@sim, PlotLegends -> {"Avg of Simulation"}, PlotStyle -> {Red, Dashed}]]]  1 answered Oct 19 '17 at 21:35 OkkesDulgerci 6,06711 gold badge1111 silver badges2121 bronze badges Here is improved version exact[\[Lambda]_, \[Mu]_, initialPopulation_] := Module[{}, DSolveValue[{p'[t] == (\[Lambda] - \[Mu]) p[t], p[0] == initialPopulation}, p[t], t]]; birthDeath[\[Lambda]_, \[Mu]_, initialPopulation_, numOfReaction_] := NestList[(\[CapitalDelta]t1 = RandomVariate[ExponentialDistribution[\[Lambda] #[[2]]]]; \[CapitalDelta]t2 = RandomVariate[ExponentialDistribution[\[Mu] #[[2]]]]; {Min[\[CapitalDelta]t1, \[CapitalDelta]t2], If[\[CapitalDelta]t1 < \[CapitalDelta]t2, #[[2]] + 1, #[[2]] - 1]}) &, {0, initialPopulation}, numOfReaction] With[{\[Lambda] = 3, \[Mu] = 1, initialPopulation = 10, numOfReaction = 200, numOfsim = 10}, sim = Table[ birthDeath[\[Lambda], \[Mu], initialPopulation, numOfReaction], numOfsim]; sim[[All, All, 1]] = Accumulate /@ sim[[All, All, 1]]; Show[ListStepPlot[sim, PlotLegends -> {"Simulation"}, PlotStyle -> Directive[AbsoluteThickness[0.2]], Frame -> True, PlotTheme -> "Detailed", FrameLabel -> {"Time", "Population"}, ImageSize -> Large], Plot[Evaluate@exact[\[Lambda], \[Mu], initialPopulation], {t, 0, Max@sim[[All, -1, 1]]}, PlotStyle -> {Black, Thick}, PlotLegends -> {"ODE"}], ListLinePlot[Mean@sim, PlotLegends -> {"Avg of Simulation"}, PlotStyle -> {Red, Dashed}]]]