Abstarct: Abstract In this thesis, we investigate the numerical solutions and behavioral analysis of the Van der Pol oscillator with a variable-order Caputo fractional derivative. First, we model the fractional variable-order Van der Pol equation and then discretize this equation using Lagrange functions and Gauss-Lobatto collocation points. Through this process, we get an algebraic system, which, when solved, yields the variables and their fractional derivatives simultaneously. We also propose an approach for calculating the fractional derivative matrix of the variable-order Caputo derivative. Additionally, we present an algorithm for solving the equation over large time intervals. By computing the residual error of our proposed method, we evaluate its high accuracy and efficiency. Next, we introduce a variable-order fractional Van der Pol equation expressed using the Caputo-Fabrizio derivative. We implement our proposed spectral-collocation method considering the properties of the Caputo-Fabrizio derivative operator and approximate the solutions and their fractional derivatives. Furthermore, the thesis explores new aspects of the variable-order fractional Van der Pol equation with time delay. First, an extended form of this equation is introduced. Then, using a spectral-collocation method baxsed on Legendre polynomials, we discretize the equation, resulting in a system of algebraic equations. Solving this system provides approximate solutions and their fractional derivatives. We also introduce a new form of coupled fractional-variable-order Van der Pol equations. A spectral-collocation approach is used to discretize the equation and compute approximate solutions and their fractional derivatives. This method is implemented over large time intervals, and the residual error calculation confirms its accuracy and correctness. At the end of each chapter, approximate solutions, phase portraits, and residual errors for different cases are presented.