I am specifically trying to find a source of heat, with heat sensors giving readings from 4 different points eg below:

The diagram above represents a 4 x 4 grid where the numbers are heat sensors giving temperature readings and each point is 100 metres apart, if we assume that normal air temp is 20, how can I find the exact point on this grid (fractional distances between points are fine) where the source of heat is?
You will need to specify a lot more information about the problem.
Are your sensors reporting the heat at steady state? At some finite time after heat started flowing from the point? Is that time known?
What is a heat source? Is it a single infinitesimal point? Several discrete points? A curve/area? Is it a source in that it is keeping the temperature constant there even as the heat diffuses away? Or is it injecting heat at a constant rate?
What is the role of the "normal air temp"? Is it just the initial temperature?
What are your boundary conditions?
The grid is made out of air? Are you assuming only diffusive heat flow (the simplest possible setting)? Or are you also taking into account convection? If so, what are you assuming about the Reynolds number of the air, the compressibility, etc? (NB in any case involving convection, finding a solution is likely utterly hopeless without sophisticated numerical simulation).
This problem is, for most formulations, a very challenging inverse problem. There are a few cases when it can be easily solved; here's one of them, adhering as much as possible to the conditions clarified below.
Let's assume that $u(x,y,t)$ is the temperature of a point $(x,y)$ on the plane at time $t$, and that there is a single heat source in the shape of a disk of radius $\epsilon$ centered at $(x_s, y_s)$ (making the heat source a disk rather than a point avoids significant complications). Let's assume that the heat at this heat source is held fixed at $u_s$, and that the heat, at a distance $d$ away from $(x_s, y_s)$, is kept constant by the ambient air at a temperature $u_d$. We will further assume that heat flow is diffusive only, with no convection or radiation, so that the temperature at any point on the annulus $\Omega$ with inner radius $\epsilon$ and outer radius $d$ is given by the heat equation
\begin{align*} \frac{du}{dt}(x,y) &= k\Delta u(x,y), \quad \epsilon \leq \|(x,y)-(x_s,y_s)\| \leq d\\ u(x, y) &= u_s, \quad \|(x,y)-(x_s,y_s)\|=\epsilon\\ u(x, y) &= u_d, \quad \|(x,y)-(x_s,y_s)\|=d, \end{align*} where $k$ is a conductivity coefficient. We are interested in the temperature at steady state $t\to \infty$, where $k$ doesn't matter and neither does the initial temperature $u(x,y,0)$. One can show that the steady state temperature is of the form
$$u(x,y) = A \log\|(x,y)-(x_s,y_s)\| + B$$ for some unknown coefficients $A$ and $B$. Plugging in the boundary conditions yields
$$u(x,y) = \frac{u_d-u_s}{\log d-\log \epsilon} \log\frac{\|(x,y)-(x_s,y_s)\|}{\epsilon} + u_s.$$
Now to the fitting. You have a list of sensor data $x_i, y_i, u_i$ and want to recover as best as possible the unknown $u_s, x_s, y_s$. This can be done by solving the nonlinear least squares problem $$\min_{x_s, y_s, u_s} \sum_i \|u(x_i,y_i)-u_i\|^2.$$