Take for example $$\int 2x \cos(x^2)dx$$ Its easy to see that this is the result of chain rule. If we take $u=x^2$ then $dx = du/2x$ and then we get $$\int \cos(u)du$$ and it is simple from there.
But this method requires algebraically manipulating the differentials, which I am aware is rigorously proven in nonstandard analysis, but how else would we arrive to this simpler form? That is when we are treating $dx$ just as notation.
My main confusion is how to justify the change of variable of integration when we are not just substituting something in place of $dx$
There really are no differentials involved here! I will use a slightly non standard notation: If $f:[a,b]\to\mathbb{R}$ is integrable i simply write $\int_a^b f$ for the Riemann integral of $f$ over $[a,b]$ instead of the more usual "$\int_a^b f(x) dx$". In what follows there will be not a single differential.
The chain rule (for the Riemann integral) says (i.e. here):
Example: Say we want to calculate the following integral $I$ defined by: $$ I:= \int_0^{\sqrt{\pi/2}} x \mapsto 2x \cos (x^2).$$
Define $g:[0,\sqrt{\pi/2}]\to [0,\pi/2]$ by $g(x) := x^2$. Then $g$ is continuously differentiable with derivative $g^\prime :[0,\sqrt{\pi/2}]\to \mathbb{R}$ given by $g^\prime(x) = 2 x$. Now define $f: [0, \pi/2]\to \mathbb{R}$ by $f(y) := \cos(y)$. Clearly $f$ is continuous.
By definition $$I = \int_0^{\sqrt{\pi/2}} g^\prime \cdot f \circ g$$ and so using the above theorem and the well known anti-derivative of $y \mapsto \cos(y)$: $$ I = \int_{g(0)}^{g(\sqrt{\pi/2})}f = \int_{0}^{\pi /2} f = \int_{0}^{\pi /2} y \mapsto \cos (y) = \sin (\pi/2) - \sin (0) = \sin (\pi/2) =1. $$