Given a function $u$ defined on $\mathbb{R}^n$, its fractional $p$-Laplacian is given by
$$
(-\Delta)_p^s u(\mathbf{x})
= C_{n,s,p}\,\mathrm{P.V.}\!\int_{\mathbb{R}^n}
\frac{|u(\mathbf{x})-u(\mathbf{y})|^{p-2}\bigl(u(\mathbf{x})-u(\mathbf{y})\bigr)}
{\|\mathbf{x}-\mathbf{y}\|^{n+sp}}d\mathbf{y},
\qquad \mathbf{x}\in\mathbb{R}^n,
$$
where the integral is understood in the principal value sense, $p \in (1,\infty)$, $s \in (0,1)$, and $C_{n,s,p}$ is a normalization constant depending on $n$, $s$, and $p$.

An equivalent formulation is given by the nonlinear Balakrishnan formula:
$$
(-\Delta)_p^s u(\mathbf{x})
= \tilde C_{n,s,p}
\int_0^\infty
\Delta (t-\Delta)^{-1}
\bigl[ \Phi_p(u(\mathbf{x}) - u(\cdot)) \bigr](\mathbf{x})
\frac{dt}{t^{1-sp/2}},
$$
where $\tilde C_{n,s,p}$ is another normalization constant, and $\Phi_p(t) = |t|^{p-2}t$. In the linear case $p = 2$, this reduces to the classical fractional Laplacian.

In this talk, we present a matrix-based spectral method for the numerical approximation of both the fractional Laplacian and the fractional $p$-Laplacian for functions defined on $\mathbb{R}^n$. Our approach builds on the Balakrishnan representation, where we discretize the second-order derivatives in $\Delta$ using spectrally accurate differentiation matrices. A key advantage is that these matrices can be diagonalized in a well-conditioned manner, enabling a stable and robust numerical scheme that naturally extends to arbitrary spatial dimensions $n$.

The method avoids multidimensional quadrature, domain truncation, and variational formulations, making it both computationally efficient and conceptually straightforward.