Composition of mappings not working as expected

Your functions take two arguments as input. Yet they return one argument as output. So the output of f cannot serve as the input of g and vice versa. You can't compose the functions with until the outputs of one function can serve as the input for the other function. Reduce your arguments in f, g by placing them in a list.

This appears to work:

f[{x_, y_}] := {E^x + y, Sin[2 x]}
g[{x_, y_}] := {2 x + Cos[y], E^(x + y)}
h[{x_, y_}] := g[f[{x, y}]]

Example

g[f[{2, 3}]]
f[g[{2, 3}]]

When f and g have been defined this way :

f[x_, y_] := {E^x + y, Sin[2 x]}
g[x_, y_] := {2 x + Cos[y], E^(x + y)}

the problem is that g and h being two-argument functions, then g is called as one-argument (i.e. a list returned by f) function h[x_, y_] := g[f[x, y]]. Instead of redefining f and g you can try another definitions of h as a composition of f and g, e.g. h1 using Sequence (and Apply) or h2 using Apply only :

h1[x_, y_] := g[ Sequence @@ f[x, y]]
h2[x_, y_] := g @@ f[x, y]

they both work well :

h1[x, y] == h2[x,y]
h1[x, y]

True

A little late to the party, but here's a function I wrote to compose arbitrary $\mathbb{R}^M \to \mathbb{R}^P$ functions with $\mathbb{R}^P \to \mathbb{R}^N$ functions:

Clear[multiDimComposition]
multiDimComposition[flst__]:=
 With[{fcns = Reverse@List[flst]},
  Fold[#2[ Sequence @@ #1 ]&, First[fcns][##], Rest[fcns]]&
 ]

Here's a few examples

Clear[f, g]
f[x_, y_] := Sin[2 \[Pi] x y^2]
g[s_] := {s, s^3}
multiDimComposition[f, g][s]
(*  Sin[2 Pi s^7] *)

The reason I created it was to specify arbitrary paths in multidimensional functions, as follows

GraphicsRow[{Show[
   DensityPlot[f[x, y], {x, -1, 1}, {y, -1, 1}], 
   ParametricPlot[g[s], {s, -1, 1}, PlotStyle -> Black]
   ], Plot[multiDimComposition[f, g][s], {s, -1, 1}]}]

Here's the transformation to spherical coordinates:

Clear[f, g]
f[x_, y_, z_] := Exp[Sqrt[x^2 + y^2 + z^2]]
g[r_, t_, f_] := {r Sin[t] Cos[f], r Sin[t] Sin[f], r Cos[t]}
multiDimComposition[f, g][Rho, Theta, Phi] // 
   Simplify[#, Rho > 0] &
(* E^Rho *)

There is a flaw, though, as written the functions must have downvalues, or it won't work:

multiDimComposition[f, {s, s^2}][s]
(* f[s] *)