Kinetic effects of bulk and fast particles in tokamaks can significantly impact the stability predictions of Magnetohydrodynamic (MHD) modes. Hybrid kinetic-MHD models are typically used to estimate these effects. Many such existing codes compute the kinetic effects by solving a reduced drift-kinetic or gyrokinetic equation semi-analytically, which often involves approximations limiting the type of orbits that can be accurately reproduced. Others solve particle orbits as a time evolution problem, requiring less assumptions but strongly increasing numerical cost. In this thesis, we develop a new spectral hybrid kinetic-MHD code, VENUS-KMHD, solving the MHD and drift-kinetic equations self-consistently. Using the Case-Van Kampen approach rather than the Landau approach, the linearised equations form a standard generalised eigenvalue problem solvable by regular iterative solvers. The considered model involves no approximation on the guiding centre trajectories of particles, thus allowing for full finite orbit width effects. The model also allows for strong flows as well as diamagnetic and perturbed parallel electric effects, whose implementation is left for future work. Our numerical method will be able to capture all these contributions while retaining the structure of the discretised problem. A first version of the new code is presented, in particular, core numerical methods and implementation details. We discuss numerical performance and challenges, key numerical parameters and the parallel implementation of the code. Detailed benchmarks of VENUS-KMHD against analytical derivations in cylindrical geometry are shown for 1/1 internal kink modes. The numerical solution of the drift-kinetic equation and resulting kinetic perturbed pressures match analytical results. Governing equations for the kinetic-MHD problem are derived in both low and high frequency limits and used to successfully benchmark the code. In the low frequency case, the ideal destabilising adiabatic effect due to hot ions is exactly cancelled by their stabilising kinetic effects. In the high frequency limit, only a weak stabilising kinetic effect is observed, requiring very accurate analytical kinetic corrections to properly verify the numerical results. A general case in the cylinder outside the high and low frequency limits is considered. Examples of stabilising and destabilising kinetic effects are presented. In a cylinder, they are often stabilising and weak but can become destabilising with an off-axis peaked kinetic pressure profile. The VENUS-KMHD numerical approach straightforwardly handles the wave-particle resonances. The ability of the code to correctly reproduce several subdominant modes is illustrated. The mode structure of the dominant and first subdominant modes is successfully benchmarked in the large frequency limit against the governing equation and an analytical result. This thesis is closed with preliminary results in toroidal geometry. The perturbed distribution function is studied in velocity space, discussing trapped particle prevalence across various hot particle temperatures. A first comparison with analytical results is presented in the high frequency (Kruskal-Oberman) limit, showing that VENUS-KMHD correctly reproduces this limit, including kinetic corrections to the inertia. Crucially, the code recovers expected kinetic corrections to the inertia in a torus, which is a significant challenge for Landau-based eigenvalue codes.