  Radial basis functions are typically used when discretization sche-mes
require inhomogeneous node distributions. While spawning from a desire to
interpolate functions on a random set of nodes, they have found successful
applications in solving many types of differential equations. However, the
weights of the interpolated solution, used in the linear superposition of basis
functions to interpolate the solution, and the actual value of the solution are
completely different. In fact, these weights mix the value of the solution with
the geometrical location of the nodes used to discretize the equation. In this
paper, we used nodal radial basis functions, which are interpolants of the
impulse function at each node inside the domain. This transformation allows to
solve a linear hyperbolic partial differential equation using series expansion
rather than the explicit computation of a matrix inverse. This transformation
effectively yields an implicit solver which only requires the multiplication of
vectors with matrices. Because the solver requires neither matrix inverse nor
matrix-matrix products, this approach is numerically more stable and reduces
the error by at least two orders of magnitude, compared to other solvers using
radial basis functions directly. Further, boundary conditions are integrated
directly inside the solver, at no extra cost. The method is naturally
conservative, keeping the error virtually constant throughout the computation.
