Question: Generalising to arbitrary curvilinear orthogonal coordinate systems

with(ArrayTools):
Operator:= proc(F::Array,X::Array) local f_11, f_12, f_21, f_22, f, P, S;

f_11:= diff(F[1],X[1]);

f_12:= diff(F[1],X[2]);

f_21:= diff(F[2],X[1]);

f_22:= diff(F[2],X[2]);

P:=Matrix(2,2,[f_11,f_21,f_12,f_22]);

S:=Matrix(2,1,[diff(G(x,y),x), diff(G(x,y),y)]);

if Size(P,2) = Size(S,1)

then return P.S;

else return

false;

end if;

end proc;

My procedure so far (I think) computes the deformation gradient and returns the transformed gradient. I am now tasked with generalisng this to arbitrary curvilinear orthogonal coordinate systems. I am unsure of how to do this in maple and what packages are available to do this. Would I simply be able to input a new line that would allow these types of coordinates, or would i have to involve a transformation that turns the coordinates into Cartesian? I also ideally want to make it work for N-dimensions, but I am unsure how to do an array with n variables