I am not sure that I fully understand the question, but in SageMath you may do the following:
sage: M = Manifold(2, 'M') # the 2-dimensional plane as a smooth manifold
sage: X.<s,t> = M.chart() # parabolic coordinates
sage: g = M.riemannian_metric('g') # the standard Euclidean metric in the plane
sage: g[0,0] = s^2 + t^2 # the metric coefficient g_{ss}
sage: g[1,1] = s^2 + t^2 # the metric coefficient g_{tt}
sage: g.display()
g = (s^2 + t^2) ds*ds + (s^2 + t^2) dt*dt
sage: F = M.scalar_field(function('f')(s,t)) # a generic smooth function M --> R
sage: nabla = g.connection() # the Levi-Civita connection associated with g
sage: lapF = nabla(nabla(F).up(g)).trace() # the Laplacian of F
sage: lapF.display()
M --> R
(s, t) |--> (d^2(f)/ds^2 + d^2(f)/dt^2)/(s^2 + t^2)
sage: lapF.expr() # the coordinate expression of Lap(F) as a symbolic expression
(diff(f(s, t), s, s) + diff(f(s, t), t, t))/(s^2 + t^2)
sage: lapF.coord_function() # same thing but with abriged notations
(d^2(f)/ds^2 + d^2(f)/dt^2)/(s^2 + t^2)