Got it, @dokken . Find below a working example which shows the same.

```
# Imports
from dolfinx import default_scalar_type
import dolfinx
from dolfinx.fem import (Constant, Function, FunctionSpace, assemble_scalar, form, locate_dofs_topological)
from dolfinx.fem.petsc import (LinearProblem, assemble_vector, assemble_matrix, create_vector,
apply_lifting, set_bc)
from dolfinx.io import XDMFFile
from dolfinx.mesh import create_rectangle, locate_entities, locate_entities_boundary, meshtags, CellType
from dolfinx.plot import vtk_mesh
from mpi4py import MPI
from ufl import (FacetNormal, Measure, SpatialCoordinate, TestFunction, TrialFunction, div, dot, dx, grad,
inner, lhs, rhs)
from petsc4py import PETSc
from mpi4py import MPI
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pyvista
# Mesh
nx, ny = 50, 50
mesh = create_rectangle(MPI.COMM_WORLD, [np.array([0,0]),np.array([5,5])], [nx, ny], CellType.triangle)
V = FunctionSpace(mesh, ("Lagrange",1))
# Temporal definitions
t = 0 # Start time
t_end = 1 # End time
num_steps = 1 # Number of time steps
dt = (t_end-t)/num_steps # Time step size
# Initialization
x = SpatialCoordinate(mesh)
T_n = Function(V)
T_n.interpolate(lambda x: np.full((x.shape[1],), 298))
T_h = Function(V)
T_h.interpolate(lambda x: np.full((x.shape[1],), 298))
T_inf = Constant(mesh,default_scalar_type(298)) # Ambient temperature
Q_gen = Constant(mesh,default_scalar_type(1000)) # Heat generation inside CV
q_tabs = Constant(mesh,default_scalar_type(500)) # Heat flux to tabs
rho = Constant(mesh,default_scalar_type(25)) # Density
C_p = Constant(mesh,default_scalar_type(90)) # Specific heat capacity of the cell
k = Constant(mesh,default_scalar_type(25)) # Thermal conductivity
h = Constant(mesh,default_scalar_type(10)) # Convection coefficient
# Creating boundaries
boundaries = [(1, lambda x: np.isclose(x[0],0)),
(2, lambda x: np.isclose(x[1],5)),
(3, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],4),np.less(x[1],5))),
(4, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],3),np.less(x[1],4))),
(5, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],2),np.less(x[1],3))),
(6, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],1),np.less(x[1],2))),
(7, lambda x: np.logical_and(np.isclose(x[0],5),np.greater_equal(x[1],0),np.less(x[1],1))),
(8, lambda x: np.isclose(x[1],0))]
facet_indices, facet_markers = [], []
fdim = mesh.topology.dim - 1
for (marker, locator) in boundaries:
facets = locate_entities(mesh, fdim, locator)
facet_indices.append(facets)
facet_markers.append(np.full_like(facets, marker))
facet_indices = np.hstack(facet_indices).astype(np.int32)
facet_markers = np.hstack(facet_markers).astype(np.int32)
sorted_facets = np.argsort(facet_indices)
facet_tag = meshtags(mesh, fdim, facet_indices[sorted_facets], facet_markers[sorted_facets])
mesh.topology.create_connectivity(mesh.topology.dim-1, mesh.topology.dim)
ds = Measure("ds", domain=mesh, subdomain_data=facet_tag)
# Variational formulation
u, v = TrialFunction(V), TestFunction(V)
a = form(rho * C_p * u * v * dx + dt * inner(k * grad(u), grad(v)) * dx +
dt * h * u * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))
L = form((rho * C_p * T_n + dt * Q_gen) * v * dx - dt * q_tabs * v * (ds(4) + ds(6)) +
dt * h * T_inf * v * (ds(1) + ds(2) + ds(3) + ds(5) + ds(7) + ds(8)))
# Solving
A = assemble_matrix(a, bcs=[])
A.assemble()
b = create_vector(L)
solver = PETSc.KSP().create(mesh.comm)
solver.setOperators(A)
solver.setType(PETSc.KSP.Type.PREONLY)
solver.getPC().setType(PETSc.PC.Type.LU)
for i in range(num_steps):
t += dt
# Update the right hand side reusing the initial vector
with b.localForm() as loc_b:
loc_b.set(0)
assemble_vector(b, L)
# Apply Dirichlet boundary condition to the vector
apply_lifting(b, [a], [[]])
b.ghostUpdate(addv=PETSc.InsertMode.ADD_VALUES, mode=PETSc.ScatterMode.REVERSE)
set_bc(b, [])
# Solve linear problem
solver.solve(b, T_h.vector)
T_h.x.scatter_forward()
# Update solution at previous time step (u_n)
T_n.x.array[:] = T_h.x.array
# Plotting the solution
pyvista.start_xvfb()
pyvista_cells, cell_types, geometry = vtk_mesh(V)
grid = pyvista.UnstructuredGrid(pyvista_cells, cell_types, geometry)
grid.point_data["u"] = T_n.x.array
grid.set_active_scalars("u")
plotter = pyvista.Plotter()
plotter.add_text("uh", position="upper_edge", font_size=14, color="black")
plotter.add_mesh(grid, show_edges=True)
plotter.view_xy()
if not pyvista.OFF_SCREEN:
plotter.show()
else:
figure = plotter.screenshot("robin_neumann_dirichlet.png")
```

The solution I got is: