  Spectral and spectral element methods using Galerkin type formulations are
efficient for solving linear fractional PDEs (FPDEs) of constant order but are
not efficient in solving nonlinear FPDEs and cannot handle FPDEs with
variable-order. In this paper, we present a multi-domain spectral collocation
method that addresses these limitations. We consider FPDEs in the
Riemann-Liouville sense, and employ Jacobi Lagrangian interpolants to represent
the solution in each element. We provide variable-order differentiation
formulas, which can be computed efficiently for the multi-domain discretization
taking into account the nonlocal interactions. We enforce the interface
continuity conditions by matching the solution values at the element boundaries
via the Lagrangian interpolants, and in addition we minimize the jump in
(integer) fluxes using a penalty method. We analyze numerically the effect of
the penalty parameter on the condition number of the global differentiation
matrix and on the stability and convergence of the penalty collocation scheme.
We demonstrate the effectiveness of the new method for the fractional Helmholtz
equation of constant and variable-order using $h-p$ refinement for different
values of the penalty parameter. We also solve the fractional Burgers equation
with constant and variable-order and compare with solutions obtained with a
single domain spectral collocation method.
