Changing the domain of a 3D Finite Difference code from cube to sphere

I have an explicit FD (Finite Difference) code for diffusion/heat on a PDE in a cuboid domain, and it works fine. I would like to update the discretized equations and change the code so as to solve the problem on a sphere. From some preliminary searches, I can't seem to find any examples of FD on spheres using rectangular coordinates.

I have seen some in spherical coordinates, but that may require too many changes to the code (if this is the only way, please advice on how best to proceed)

$$\frac{\partial U}{\partial t} = \alpha \nabla^2U\quad in\enspace\mathbb{R}^3$$ $$U(x,y,z,0)=T$$ where T is a constant, and is also the boundary condition applied on all faces of the cube / the surface area of the sphere.

Question: Is it possible to adapt my code for a cube, to work on a sphere, and remain in rectangular coordinates? If so, how do I approach this?

Edit: (more context) This is a homogenous heat equation, where the boundary condition is a constant temperature (say T=100 degrees) applied on the entire surface area of the sphere (for the cube it was also applied on all 6 faces). The initial condition is the same heat source on the surface, while the entire insides of the sphere was initially at zero.

That is certainly possible. You can reformulate your problem to use $$\Delta + V$$ instead of $$\Delta$$ as an operator, where $$V=V(x,y,z)=V(\vec{r})$$ is a potential defined in a way suitable to tackle your problem. For instance, if I were looking for eigenmodes of a disc or a sphere by solving the Helmholtz equation using finite differences I'd define the potential as
$$V(\vec{r}) = \begin{cases}0 & \vec{r}< r_0\\ V_0 & else\end{cases}$$
with $$V_0\rightarrow\infty$$ being essentially a (very) large number. This should, in the case of a Helmholtz problem, force the resulting solution to be ~0 outside of the sphere with radius $$r_0$$. For a heat equation you can do something similar (depending on the boundary condition).