I am having trouble understanding the similarity solution method for solving partial differential equations. I have been able to replicate the highly spoon-fed example, but all of the other, more general explanations online introduce too many stipulated variables and contain too much unmotivated computation for me to follow. Say I want to solve $u_t(x,\ t) = 3u_{xx}(x,\ t)$. I can go back and add conditions in when I need them.
I think I understand the first step, albeit in a more hand-wavy manner than is generally used. The goal is to find constants $a$, $b$, and $c$ such that $u(x,\ t)$ being a solution to the PDE implies that $e^c u(e^a x,\ e^b t)$ is also a solution to the PDE. This is the scaling transformation $x \rightarrow e^a x, t \rightarrow e^b t, u \rightarrow e^c u$. I plug these into the "partial differentials" of the PDE and pull out the constants, as one does when nondimensionalizing a system. This gives me that $e^{c - b}u_{t} = 3e^{c - 2a}u_{xx}$. Thus, in order for this to be arithmetically consistent with the original PDE, it must be the case that $2a = b$, while $c$ is arbitrary. So now we have that $u(x,\ t)$ and $e^c u(e^a x,\ e^{2a} t)$ are similar solutions.
I am lost on the logic for progressing past this point. I don't think linking an online explanation will be much help, as I have already looked into many of them. What is the rationale for the remainder of the solution method from this point, and how can I continue the computation in an intuitive way?
Here is one way. I'll write $y=px, \tau = qt$ rather than using the exponential notation, as you say like changing units in physics, which to me is the main motivation for this method. I'm sure there are more sophisticated points of view.
Now suppose $u$ solves the heat equation $u_t = u_{xx}$ and $$ u(x,t) = rv(px,qt). $$ Then as you know, $$ u_t-u_{xx} = rqv_\tau-rp^2 v_{yy} $$ shows that $v$ is also a solution whenever $q = p^2$.
A new idea: suppose $v$ is not just some new solution, but is actually $u$ itself, $$ u(x,t) = u(px,p^2 t) $$ for all $x,t,p$. [This is of course not likely! But we check it.] We can then eliminate the $t$ variable by setting $p = t^{-1/2}$: $$ u(x,t) = u(t^{-1/2}x,1). $$ To see whether there really is a solution that depends only on $t^{-1/2}x$, we try $$ u(x,t) = f(t^{-1/2}x) $$ in the heat equation. Denote $s=t^{-1/2}x$. We get $$ u_t-u_{xx} = -\tfrac{1}{2}t^{-3/2}xf'(s)-t^{-1}f''(s) = -t^{-1}\big(\tfrac{1}{2}sf'(s)+f''(s)\big). $$ So this works if we have $$ f''+\tfrac{1}{2}sf' = 0. $$ Solving that you get $$ f(s) = c_1+c_2\int_{-\infty}^s e^{-w^2/4}dw $$ or $$ u(x,t) = c_1+c_2\int_{-\infty}^{t^{-1/2}x} e^{-w^2/4}dw. $$ To clean that up a little, note that whenever $u$ is a solution to the heat equation, then so is any multiple of $u_x$. That gives a famous solution $$ t^{-1/2}e^{-x^2/4t}. $$