gridtransform := proc( u, v, xr, yr, n, ax)

   local i, j, x, y, p, q, r, s, dx, dy, L1, L2, k, a, b, c, d:
   
   a := op(1,xr); b := op(2,xr);
   c := op(1,yr); d := op(2,yr);
   p := array(0..n):  q := array(0..n):
   r := array(0..n):  s := array(0..n):
   x := array(0..n+1):  y := array(0..n+1):
   x[0] := a:  y[0] := c:
   dx := evalf( (b - a) /n):  dy := evalf( (d - c) / n):
   
   for i from 0 to n do
      r[i] := plot([ x, y[i], x = a..b]):
      s[i] := plot([x[i] ,y, y = c..d]):
      p[i] := plot( [ u(x[i], y), v(x[i], y), y = c..d]):
      q[i] := plot( [ u(x, y[i]), v(x, y[i]), x = a..b]):
      x[i + 1] := x[i] + dx:
      y[i + 1] := y[i] + dy:
   od:

   L1 := display({ r[k] $ k = 0..n, s[k] $ k = 0..n}, 
         scaling =constrained, axes = ax): 
   L2 := display( { p[j] $ j = 0..n, q[j] $ j = 0..n}, 
         scaling = constrained, axes = ax ): 
   display([ L1, L2 ], insequence = true);
   
end: