💼 Case: consider the potential $P(\vec r)$ at $\vec r$ due to a charge $q$ at $\vec r'$
$$ \begin{aligned} \vec R&=\vec r-\vec r' \\ &=(x-x')\,\hat x+(y-y')\,\hat y \\ &\quad+(z-z')\,\hat z \end{aligned} $$
$$ \hat R=\frac{\vec r- \vec r'}{|\vec r- \vec r'|} =\frac{\vec R}{|\vec R|}=\frac{\vec R}{R} $$
https://www.desmos.com/calculator/hgwq2ryly2
$$ R=\underbrace{\left [ (x-x')^2 +(y-y')^2+(z-z')^2 \right ]^\frac 12 }{\text{cartesian}}=\underbrace{\left [r^2+r'^2-2rr'\cos(\alpha)\right ]^\frac 12}{\text{angle form}} $$
$$ \boxed{V(\vec r)=\frac{q}{4\pi \epsilon_0} \frac{1}{R}} $$
Which we can use to find the electric field
$$ \begin{aligned} \vec E(\vec r)&=-\vec \nabla V(\vec r)=-\frac{q}{4\pi \epsilon_0} \vec \nabla \left ( \frac{1}{R} \right )=\frac{q}{4\pi\epsilon_0}\frac{\hat R}{R^2} \end{aligned} $$
🗒️ Note: $x$ and $x'$ are independent so $\frac{\partial x'}{\partial x}=\frac{\partial y'}{\partial y}=\frac{\partial z'}{\partial z}=0$
For a set of point charges the potential at each point will be the sum of the potential of each charge
$$ \begin{aligned} V(\vec r)&=\frac{1}{4\pi \epsilon_0} \sum_i \frac{q_i}{R_i} \qquad \text{where } \; \vec R_i = \vec r- \vec r_i'\\ \vec E(\vec r) &= \frac{1}{4\pi \epsilon_0~} \sum_i q_i \frac{\hat R_i}{R_i^2} \end{aligned} $$
Line charge
$$ \text d q'=\lambda (\vec r')\,\text d l' $$
Surface charge
$$ \text d q'=\sigma(\vec r')\,\text da' $$
Volume charge
$$ \text d q' = \rho (\vec r') \, \text d \tau' $$
$$ \begin{aligned} V(\vec r)&= \frac{1}{4\pi \epsilon_0} \int \frac{ \lambda (\vec r') \, \text d l'}{R} \\ \vec E ( \vec r)&= \frac{1}{4\pi \epsilon_0} \int \frac{ \hat R}{R^2}\lambda (\vec r') \, \text d l' \end{aligned} $$
$$ \begin{aligned} V(\vec r)&= \frac{1}{4\pi \epsilon_0} \int_S \frac{ \sigma(\vec r') \, \text d a'}{R} \\ \vec E ( \vec r)&= \frac{1}{4\pi \epsilon_0} \int_S \frac{ \hat R}{R^2}\sigma (\vec r') \, \text d a' \end{aligned} $$
$$ \begin{aligned} V(\vec r)&= \frac{1}{4\pi \epsilon_0} \int_V \frac{ \rho(\vec r') \, \text d \tau'}{R} \\ \vec E ( \vec r)&= \frac{1}{4\pi \epsilon_0} \int_V \frac{ \hat R}{R^2}\rho(\vec r') \, \text d \tau' \end{aligned} $$
Starting from the electric field
$$ \vec E=-\vec \nabla V $$
we can write
$$ \begin{aligned} &\;\vec \nabla \times \vec E=-\vec \nabla \times (\vec \nabla V)\\ &\boxed{\vec \nabla \times \vec E=0} \end{aligned} $$
Flux $\Phi$ of $\vec E$ through an arbitrary closed surface $S$ that encloses a point charge
Define $\beta$ to be the normal angle between $\theta$ and the surface that way $\text da$ can represent any shape
$$ \text da=\frac{R\sin(\theta)\,\text d \theta \text d \phi}{\cos(\beta)} $$
Now we calculate the flux
$$ \begin{aligned} &\text d \Phi=\vec E \cdot \text d \vec a =E\,\text da\, \cos(\beta) \\ &\;\quad =\frac{q}{4\pi \epsilon_0}\sin(\theta) \, \text d \theta \text d \phi \\ &\;\; \Phi=\frac{q}{4\pi\epsilon_0}\int_0^\pi \sin(\theta)\,\text d \theta\int_0^{2\pi} \,\text d\phi \\ &\;\quad=\frac{q}{4\pi \epsilon_0}[1-(-1)][2\pi-0] \\ &\;\boxed{\Phi=\frac{q}{\epsilon_0}} \end{aligned} $$
Using the principle of superposition, the result for a surface enclosing multiple point charges
$$ \Phi=\frac{1}{\epsilon_0} \sum_i q_i $$