3D equation. You can write a PDE for conduction in solids with continuous properties. This equation comes from the energy balance equation and reads
$\dfrac{\partial e}{\partial t} = - \nabla \cdot \mathbf{q}$
being $e$ the energy that (under some assumptions) can be written as $e = \rho c T$ and $\mathbf{q}$ the heat flux vector, that can be written as $\mathbf{q} = - k \nabla T$ if Fourier's law holds. To complete the problem, you need initial and boundary conditions.
1D equation. Let's focus now on the 1D dimensional equation we can derive for elongated solid elements, like beams or rods.
- assuming constant density $\rho$ and heat capacity $c$ (that we thus take outside the time derivative),
- after integration on the section of the solid element
$\rho c A \dfrac{\partial T}{\partial t} = - \dfrac{\partial (A q)}{\partial x} - Q_{lateral}$,
being $Q_{lateral}$ the heat flux exchanged through the lateral surface of the element. Assuming that the element lateral surface is insulated, $Q_{lateral} = 0$. The 1D expression of Fourier's law reads
$q = - k \dfrac{\partial T}{\partial x}$,
and thus the differential equation reads
$\rho c A(x) \dfrac{\partial T}{\partial t}(x,t) = \dfrac{\partial}{\partial x} \left( A(x) k(T) \dfrac{\partial}{\partial x} \left( T(x,t) \right) \right)$$\quad, \quad t \in [0, T_{fin}], \quad x \in [x_1, x_2]$
with:
- the initial conditions $T(x,0) = T_0(x)$, $x \in [x_1, x_2]$
- the boundary conditions in $x_1$, $x_2$ for $t \in [0, T_{fin}]$, as an example:
- prescribed temperature (Dirichlet or essential boundary conditions): $T(x_1, t) = T_1(t)$
- prescribed heat flux (Neumann or natural boundary conditions):
$Aq(x_2, t) = A(x_2) k(T) \dfrac{\partial T}{\partial x}\bigg|_{x_2} = Aq_2(t)$.
Numerical solution. Despite this equation looks quite complex, it can be solved quite in a easy way using its weak formulation and finite element approach. We can get the weak formulation as
$\displaystyle \int_{\Omega} w(x) \rho c A \dfrac{\partial T}{\partial t} + \int_{\Omega} A k \dfrac{\partial w}{\partial x}\dfrac{\partial T}{\partial x} = w(x_2) Q_2(t)$, $\forall w(x) \in W$.
(This is just a hint about the numerical method, without explicitly referring to mathematical details, especially about the weak formulation and the test space $W$ the test function $w(x)$ belongs to, and the treatment of the essential boundary conditions)