This chapter introduces key concepts in the numerical approximation of first-order partial differential equations of hyperbolic type. We begin with the simplest method, the finite difference method, applied to the basic linear advection equation. We derive several well-known schemes for hyperbolic systems and discuss their essential properties, including Godunov’s upwind, Godunov’s centred, Lax-Friedrichs, Lax-Wendroff, Fromm, Warming-Beam, and FORCE schemes. Theoretical concepts for analysing numerical methods are introduced, such as local truncation error, order of accuracy, modified equation, linear stability, and monotonicity. We also present a numerical version of the von Neumann method for studying the linear stability of numerical methods, which is applicable to high-order methods in multiple space dimensions. A shortcut for analysing accuracy is introduced, simplifying the evaluation of a numerical method’s accuracy, even in multiple space dimensions. Godunov’s theorem is stated and proven, demonstrating that methods of high-order accuracy must be nonlinear, even when applied to linear equations. The class of conservative methods is introduced, and properties of the numerical flux are examined. Classical finite difference methods are reformulated in conservative form by identifying the appropriate numerical flux. Boundary conditions, including reflective and transmissive, are discussed. We also cover the calculation of a stable time step for explicit methods, as determined by the Courant stability condition, with an emphasis on providing reliable estimates of the maximal wave speed to ensure stability in nonlinear problems. Numerical results are presented to assess the performance of the methods. Mesh refinement, along with concepts of error, mesh-independent solutions, convergence rates, and efficiency, is discussed and demonstrated through practical examples. Finally, a list of exercises is provided at the end of the chapter.