When the graph between Gibbs energy and extent of reaction is plotted, such a curve is obtained:
I am unable to understand how we got to the concavity of this graph (which can also be seen as the monontonic nature of the graph during various segments of the reaction). And that is what I am interested in knowing. I want to see the math worked out.
The math is fairly straightforward if we can start from $$\Delta_r G = \Delta G ^\circ + RT\ln Q$$
Let's use the example of a very simple reaction $\ce{A <=> B}$ with an equilibrium constant of $K=1$ (so $\Delta G^\circ=0$) and assume that the activities of both A and B are equal to their concentrations and our starting concentration of A is 1 M.
In this simple case, $Q=\frac{[\ce{B}]}{[\ce{A}]}=\frac{\xi}{1-\xi}$, which, combined with $\Delta G^\circ =0$, gives us that $$\Delta_r G = RT\ln \frac{\xi}{1-\xi}$$
As noted on the plot you showed, $\Delta_r G$ represents the slope of the plot, so to reproduce the plot, we integrate the right-hand expression and get $$G(\xi)=RT\left(\ln(1-\xi) + \xi\ln\frac{\xi}{1-\xi}\right)$$ $$=RT\left((1-\xi)\ln(1-\xi) + \xi\ln\xi\right)$$
A plot of this function in the range $0 <\xi <1$ gives something resembling your original plot, except symmetric due to the choice of $K=1$.
ADDED: If part of your question is why $\Delta_r G = \frac{dG}{d\xi}$, we can derive that as well. The total free energy of the system is $G = G_A + G_B$. For the separate components, the molar free energies (ignoring the reference concentration term (1M) for simplicity) are $$\overline{G}_A=\mu_A^\circ + RT\ln[\ce{A}]$$ and $$\overline{G}_B=\mu_B^\circ + RT\ln[\ce{B}]$$.
Since these are molar free energies, we need to multiply each by the relevant amount. Let's assume our system is 1 L, so the concentration and the amount of A are both $1-\xi$ and for B are both $\xi$.
Substituting in, we have $$G(\xi) = G_A + G_B$$ $$= (1-\xi)\overline{G}_A + \xi\overline{G}_B$$ $$= (1-\xi)(\mu_A^\circ + RT\ln(1-\xi)) + \xi(\mu_B^\circ + RT\ln\xi)$$
Recalling that we have set $K=1$ so that $\mu_A^\circ = \mu_B^\circ$, we can do a lot of rearranging and cancelling and get the integration result from above: $$G(\xi)=RT\left((1-\xi)\ln(1-\xi) + \xi\ln\xi\right)$$
