The Method of Characteristics (MoC) is a well-known procedure used to find the numerical solution of systems of hyperbolic partial differential equations (PDEs). The main idea of the MoC is to integrate a system of ordinary differential equations (ODEs) along the characteristic curves admitted by the PDEs. In principle, this can be done by any appropriate numerical method for ODEs. In this thesis, we will examine the MoC applied to systems of hyperbolic PDEs with straight-line and crossing characteristics. So far, only first- and second-order accurate explicit MoC schemes for these types of systems have been reported. As such, the purpose of this thesis is to develop MoC schemes which are of an order greater than two. The order of the global truncation error of an MoC scheme goes hand-in-hand with the order of the ODE solver used. The MoC schemes which have already been developed use the first-order Simple Euler (SE) and second-order Modified Euler (ME) methods as the ODE solvers. The SE and ME methods belong to a larger family of numerical methods for ODEs known as the Runge--Kutta (RK) methods. First, we will attempt to develop third- and fourth-order MoC schemes by using the classical third- and fourth-order RK methods as the ODE solver. We will show that the resulting MoC schemes can be strongly unstable, meaning that the error in the numerical solution becomes unbounded rather quickly. We then turn our attention to the so-called pseudo-RK (pRK) methods for ODEs. The pRK methods are at the intersection of RK and multistep methods, and a variety of third- and fourth-order schemes can be constructed. We show that when certain pRK schemes are used in the MoC, at most a weak instability, or no instability at all, is present, and thus the resulting methods are suitable for long-time computations. Finally, we present some numerical results confirming that the MoC using third- and fourth-order pRK schemes have the desired accuracy.