In this thesis, we study fractional differential equations with Hilfer derivative operators. Solutions are approximated using a newly developed Bernstein-splines approach and subsequently applied to the Van der Pol oscillator with fractional damping. Fractional derivatives generalize differentiation to the order α ∈ (0,∞), resulting in order α differential operators. The Hilfer-derivative of order α ∈ (0,1) and type β ∈ [0,1] is one of such operators and combines two of the most commonly used operators through parametrization: the Riemmann-Liouville-derivative (β = 0) and Caputo-derivative (β=1). The choice of β results in different kernel behavior of the operator, commonly yielding singular behavior of solutions of initial value problems (IVP's) and boundary value problems (BVP's) for β ∈ [0,1). Based on existing results for IVP's, we develop a new proof for existence and uniqueness of solutions to BVP's for Hilfer-fractional derivatives. To obtain solution approximations numerically for IVP's and BVP's, a Bernstein splines solution approach is developed and implemented, providing accurate convergence results for nonlinear IVP's and BVP's in an efficient vectorized approach. Implementation difficulties for the singular behavior of solutions for β ∈ [0,1) are successfully resolved through time-domain transformation approximation techniques. Finally, the methods are applied to numerically study the behavior of the fractionally damped Van der Pol oscillator, a nonlinear equation of interest in electrical engineering and control theory. We study the approximate limit cycle, of which behavior corresponds with existing analytical results for the Caputo-derivative (β=1). Convergence of solutions can also be obtained for Hilfer type values of β ∈ [0,1), appearing to be of little influence on the long-term limit cycle. Furthermore, when forcing is applied, we observe periodic, quasiperiodic and chaotic behavior.