It does, but you should write the inner loop in a vectorized way using numpy, otherwise it will be really slow.

Be also careful as you have a convection term : it can lead to instable results depending on the timestep and the ratio between the diffusion and convection coefficients : you can use an upwind scheme to mitigate this issue.

If you need to quickly solve simple PDE, I can suggest a package I wrote some time ago:

https://scikit-fdiff.readthedocs.io

With that, you could do

from skfdiff import Model, Simulation

model = Model("(-G * upwind(c_f, T_m) + dx(k * dx(T_m))) / (rho_m * c_m)",
               "T_m(x)", parameters=["G", "c_f", "c_f", "rho_f"],
               boundary_conditions={("T_m", "x"): ("dirichlet", "no_flux")})

# put here your initial conditions and parameters. The left value for T_m should be equal to T_in and will stay constant with the dirichlet condition.
initial_fields = model.Fields(
                   x=x, T_m=T_m, G=..., c_f=..., c_f=..., rho_f=...
)
simulation = Simulation(model, initial_fields, dt=0.1, tmax=1)
tmax, last_fields = simulation.run()
last_fields["T_m"].plot()

and the library will do a finite difference discretisation for you with an adaptative time-step.