Probabilistic solutions of fractional differential and partial differential equations and their Monte Carlo simulations

The paper’s Lemma 1 correctly proves items (1)–(3) using the Laplace pair ∫_0^∞ e^{-st} g_β(x; t) dt = s^{β-1} e^{-x s^β} and standard Caputo–Laplace identities (equations (1.11) and (1.4) in the paper), and the kernel subordination g_{β α}(x; t) = ∫_0^∞ g_β(x; s) g_α(s; t) ds (equation (1.12)). These are stated and used in Lemma 1 and its proof in the manuscript itself. However, the manuscript’s Lemma 1 (4) concludes D_C^{β α} f_{β α}(t) = D_C^α(D_C^β f_β(t)) by replacing ∫_0^∞ D_C^β f_β(s) g_α(s; t) ds with D_C^α(D_C^β f_β(t)), which is not a valid step in general. In Laplace domain, these two sides have different transforms for generic f (e.g., f(t) = e^{-λ t}). This inconsistency directly contradicts the asserted equality. The paper’s own statements and identities make the fix clear: what always holds is D_C^{α β} f_{β α}(t) = ∫_0^∞ f'(x) g_{α β}(x; t) dx via subordination, but not the claimed composition D_C^α(D_C^β f_β(t)). Items (1)–(3) are fine; item (4) is wrong as stated and proved. See the definitions and identities (1.4), (1.11), (1.12) and Lemma 1 in the paper for the precise formulas cited here .