Explaining the method of characteristics

The basic idea is that we look for parametric curves along which the PDE tells us how the function changes. Suppose you have a smooth parametric curve in the $x-y$ plane: $x = X(t)$, $y = Y(t)$. Consider how a smooth function $u(x,y)$ changes as you move along the curve. The chain rule says $$\frac{du}{dt} = \frac{\partial u}{\partial x} \frac{dX}{dt} + \frac{\partial u}{\partial y} \frac{dY}{dt}$$ Now this looks rather like the left side of a first-order PDE $$a(x,y) \frac{\partial u}{\partial x} + b(x,y) \frac{\partial u}{\partial y} = c(x,y)$$ In fact, if you can find $X(t)$ and $Y(t)$ such that $\frac{dX}{dt} = a(X(t),Y(t))$ and $\frac{dY}{dt} = b(X(t),Y(t))$, this tells you that along that curve $\frac{du}{dt} = c(X(t), Y(t))$, so that if you know $u(X(0), Y(0))$ you can get $$u(X(s),Y(s)) = u(X(0), Y(0)) + \int_0^s c(X(t),Y(t))\ dt$$

Now, for any point $(x_1, y_1)$ in the plane, suppose you want to find $u(x_1, y_1)$, where $u(x,y)$ satisfies the PDE $a(x,y) \frac{\partial u}{\partial x} + b(x,y) \frac{\partial u}{\partial y} = c(x,y)$ plus some boundary condition. Then you want to find a curve $x = X(t)$, $y = Y(t)$ satisfying the system of ordinary differential equations $\frac{dX}{dt} = a(X(t),Y(t))$ and $\frac{dY}{dt} = b(X(t),Y(t))$ that passes through the given point $(x_1, y_1)$, follow that curve to some $(x_0, y_0)$ where your boundary condition tells you the value of $u$, and then you can get $u(x_1, y_1)$ by an integral as above.

This is an intuitive take which is intended to be as simple as possible.

Let's say you're measuring some system for two variables $x$ and $t$, and you discover that you can describe the system by $u_t+a(x,t)u_x=0$.

Okay, but what if $x$ and $t$ just superficially look like independent variables, but in fact $x$ is a function of $t$?

If it were so, then $$\frac{d}{dt}u(x(t),t)=u_t+x'(t)u_x$$ which means that if we can find a function $x(t)$ with the property that $x'(t)=a(x,t)$, then $u(x(t),t)$ would have the property that $$\frac{d}{dt}u(x(t),t)=u_t+a(x,t)u_x$$ which is exactly the type of function we are looking for.

So if it really is the case that our measured variable $x$ actually is a function of the underlying variable $t$, we only need to find a function $x(t)$ such that $x'(t)=a(x(t),t)$ with $x(0)=x_0$ and solve the much easier ODE $$\frac{d}{dt}u(x(t),t)=0$$ and then $u(x(t),t)$ satisfies the PDE we started with.

Usually these problems come with some boundary conditions like $u(x,0)=\phi(x)$, so we solve our easy ODE: $$u(x(t),t)=C$$

If we set $C=\phi(x_0)$ we see that $u(x(t),t)=u(x(0),0)=u(x_0,0)=\phi(x_0)$ (indeed, $u$ is constant wrt $t$)

So $$u(x(t),t)=\phi(x_0)$$

is a solution to the PDE provided we can find an $x(t)$ such that $x'(t)=a(x(t),t)$ with $x(0)=x_0$, which we can. When you find $x(t)$, you can solve it for $x_0$ and substitute so you get a solution function consisting of $x$'s and $t$'s, which you can then differentiate to check that it satisfies the PDE.