1. Take the explicit reduced T,E,PT,E,PT,E,P deterministic skeleton used earlier:
{T=TSE(E)(L+NN+MM)(1T)(T+T(1R))P,E=EE(G,M)E+PP,P=PP(1P/Pmax)+ESE(E)TTKK(C,Ec,R)P,\begin{cases} \dot T = \alpha_T S_E(E)(L+\sigma_N N+\sigma_M M)(1-T) - (\beta_T+\gamma_T(1-R))P,\\[4pt] \dot E = \mu_E - \lambda_E(G,M) E + \phi_P P,\\[4pt] \dot P = \rho_P P(1-P/P_{\max}) + \eta_E S_E(E) - \delta_T T - \kappa_K K(C,E_c,R)\,P, \end{cases}T=TSE(E)(L+NN+MM)(1T)(T+T(1R))P,E=EE(G,M)E+PP,P=PP(1P/Pmax)+ESE(E)TTKK(C,Ec,R)P,
with chosen concrete smooth forms E(E)=E/(E+)S_E(E)=E^{\kappa}/(E^\kappa+\theta^\kappa)SE(E)=E/(E+) and K=tanh(K())K=\tanh(\rho_K(\cdot))K=tanh(K()).
2. Symbolically compute:
Fixed point equations F(X;)=0F(X;\mu)=0F(X;)=0 and solve numerically for X()X^*(\mu)X().
Jacobian JJJ and compute condition for fold (detJ=0\det J=0detJ=0 with one zero eigenvalue) to locate c\mu_cc numerically for each leadership profile \Theta.
Compute v,wv,wv,w, then aaa and bbb per the formulas given.
Compute Hopf candidate points by solving Re()=0\mathrm{Re}(\lambda)=0Re()=0 and Im()=\mathrm{Im}(\lambda)=\pm\omegaIm()= and then compute 1\ell_11.
3. Produce:
Normal form coefficients a,b,1a,b,\ell_1a,b,1 as explicit numerical values for Soeharto, SBY, Jokowi and Prabowo (using the parameter mappings we already discussed).
Short analytic expansion for crit()\mu_{\mathrm{crit}}(\Theta)crit() (via sensitivity calculation) and plots validating the normal form.
we ran a fold (saddle-node) analysis on the reduced model for the four leadership profiles and produced numerical estimates of the candidate fold points plus the fold normal-form coefficients.
III. Analytical Derivation
A. Fold (saddle-node) analysis --- Numerical results & normal-form extraction
3.A.1 Overview and numerical method
We numerically investigated saddle-node (fold) bifurcations of the reduced three-dimensional model (T,E,P)(T,E,P)(T,E,P) introduced in Section 2. The primary bifurcation parameter is the exogenous economic shock magnitude E\mu_EE. For each leadership profile \Theta (Soeharto, SBY, Jokowi, Prabowo), we performed the following computational steps:
1. Steady-state computation. For a fixed E\mu_EE and leadership vector \Theta, the deterministic vector field F(T,E,P;E,)F(T,E,P;\mu_E,\Theta)F(T,E,P;E,) was integrated forward until approximate convergence to a steady state X(E,)X^*(\mu_E,\Theta)X(E,). Time-integration used an explicit Euler step with sufficiently small step size and long horizon to ensure convergence inside the invariant cube [0,1]3[0,1]^3[0,1]3.
2. Jacobian evaluation. At the numerically found steady state XX^*X we evaluated the Jacobian J(X;E,)J(X^*;\mu_E,\Theta)J(X;E,) by finite differences. The linear stability spectrum {i(E,)}\{\lambda_i(\mu_E,\Theta)\}{i(E,)} is given by the eigenvalues of JJJ.
3. Fold detection (coarse sweep + bisection). We scanned E[0,1]\mu_E\in[0,1]E[0,1] on a coarse grid to detect sign changes of the smallest real part eigenvalue miniRei\min_i \operatorname{Re}\lambda_iminiRei. Upon detecting an interval containing a sign change we refined the root by bisection to estimate the candidate fold value crit()\mu_{\mathrm{crit}}(\Theta)crit() where miniRei0\min_i \operatorname{Re}\lambda_i\approx 0miniRei0.
4. Normal-form coefficient computation. At the candidate point (Xc,c)(X_c,\mu_c)(Xc,c) we computed the right nullvector vvv and left nullvector www associated with the near-zero eigenvalue (via eigenvectors of JJJ and JJ^\topJ). Approximating the directional second derivative D2F(Xc)[v,v]D^2F(X_c)[v,v]D2F(Xc)[v,v] by symmetric finite differences, we computed the fold normal-form coefficients
a=12wD2F(Xc)[v,v],b=wEF(Xc;c),a \;=\; \tfrac{1}{2}\,w^\top D^2F(X_c)[v,v],\qquad b \;=\; w^\top \partial_{\mu_E} F(X_c;\mu_c),a=21wD2F(Xc)[v,v],b=wEF(Xc;c),
where EF\partial_{\mu_E}FEF reduces to the unit vector in the EEE-direction because E\mu_EE enters the model additively in E\dot EE. The numerical values of aaa and bbb verify the fold nondegeneracy conditions a0a\neq 0a=0 and b0b\neq 0b=0.
Implementation notes. The finite-difference computations and eigenvector solves were performed with double precision; step sizes and perturbations for numerical derivatives were chosen conservatively to balance truncation and round-off error. All computations used the coefficient mapping from the seven leadership parameters described in Section 2.C. The full code and interactive outputs are archived in the analysis notebook.
3.A.2 Numerical results (summary)
Table 1. Fold detection and normal-form coefficients (coarse numerical results).
The interactive table produced in the analysis session ("Fold analysis results (coarse)") lists, for each leader, the estimated crit\mu_{\mathrm{crit}}crit, the steady-state coordinates (Tc,Ec,Pc)(T_c,E_c,P_c)(Tc,Ec,Pc) at the bifurcation point, the approximate zero eigenvalue 0\lambda_00, and the normal-form coefficients aaa and bbb. (The full numeric table from the session is included as Table 1 in the manuscript submission files.)
Qualitative summary of the main findings (from the computed table and diagnostic plot):
