Try using the Helmholtz decomposition. Write
$${\bf F}=\nabla\times {\bf G}-\nabla\phi$$
where
$${\bf G}({\bf r})=\frac{1}{4\pi}\int\frac{\nabla'\times {\bf F}({\bf r}')}{|{\bf r}-{\bf r}'|}\,dV'$$
$$\phi({\bf r})=\frac{1}{4\pi}\int\frac{\nabla'\cdot {\bf F}({\bf r}')}{|{\bf r}-{\bf r}'|}\,dV'$$
Hopefully, the double integrals will allow exchange of integration order and you will get something reasonable.
I didn't try it myself.
EDIT:
It states that every vector field can be split into a irrotational potential contribution and a sourceless curl contribution. The formulas are the ones above (of course it's for 3D space only). It's basically just doing "there and back again": F=derivative of integral of derivative. But it does split the contributions nicely, and it sometimes simplifies the expressions.
However, unfortunately your $\vec{F}$ is not bounded (it increases quadratically with R) so the above aren't correct without an additional surface terms (see http://en.wikipedia.org/wiki/Helmholtz_decomposition). Even more, you aren't actually looking for decomposition into curl and grad, but a simple vector expression.
It's also much more simple that it looks:
$$||\vec{R}-\vec{r}||^2=\vec{R}^2+\vec{r}^2-2\vec{r}\cdot\vec{R}$$
integrating this on a closed curve, you get
$$\oint \vec{r}^2\,d\vec{r}+||R||^2\oint d\vec{r}-2\vec{R}\oint \vec{r}\,d\vec{r}$$
The first term is independent of $\vec{R}$ and is just a vector (call it $\vec{B}$). The second term is $0$ because you are on a closed curve. The last term evaluates to a dot product of $\vec{R}$ on to a vector, that is essentially a length-rescaled position of the center of gravity of the curve: $\vec{A}=l\vec{r}_{center}$.
This means that the $\times$ in your expression isn't a cross product.
The second part is straight-forward:
$$\oint_C(\vec{r}-\vec{R})^2\,d\vec{r}$$
Take the curl
$$\nabla_R\times\oint_C(\vec{r}-\vec{R})^2\,d\vec{r}$$
$$=-2\oint_C (\vec{r}-\vec{R})\times\,d\vec{r}=-4\iint\,d\vec{S}$$
where we saw that this is completely independent of $\vec{R}$ because the curve is closed ($\vec{R}$ is just a displacement of the curve).
In fact, the only constraints for the vector $\bf{n}$ are
$1.$ The vector $\bf{n}$ is a unit vector normal to the surface.
$2.$ It should have proper orientation depending on the orientation of the surrounding curve.
So, I think you may have made a mistake in the problem you solved and hence we may help you if you write it down in your question. :)
Verifying Stokes Theorem For Your Question
Your surface is enclosed by the intersection curve of the plane $x+z=1$ and the cylinder $x^2+y^2=36$ as the following figure shows.
The parametric equation of the intersection curve, the tangent vector, and the vector field are
$$\eqalign{
& {\bf{x}} = 6\cos \theta {\bf{i}} + 6\sin \theta {\bf{j}} + \left( {1 - 6\cos \theta } \right){\bf{k}} \cr
& {{d{\bf{x}}} \over {d\theta }} = - 6\sin \theta {\bf{i}} + 6\cos \theta {\bf{j}} + 6\sin \theta {\bf{k}} \cr
& F({\bf{x}}) = xy{\bf{i}} + 2z{\bf{j}} + 6y{\bf{k}} \cr} $$
and hence the line integral will be
$$\eqalign{
& I = \int\limits_C {F({\bf{x}}) \cdot {{d{\bf{x}}} \over {d\theta }}d\theta } = \int_{\theta = 0}^{2\pi } {\left( { - 6\sin \theta xy + 12\cos \theta z + 36\sin \theta y} \right)d\theta } \cr
& \,\,\, = 6\int_{\theta = 0}^{2\pi } {\left( { - 36{{\sin }^2}\theta \cos \theta + 2\cos \theta \left( {1 - 6\cos \theta } \right) + 36{{\sin }^2}\theta } \right)d\theta } \cr
& \,\,\, = 6\int_{\theta = 0}^{2\pi } {\left( { - 36{{\sin }^2}\theta \cos \theta - 12{{\cos }^2}\theta + 36{{\sin }^2}\theta + 2\cos \theta } \right)d\theta } \cr
& \,\,\, = 6\left[ { - 36\int_{\theta = 0}^{2\pi } {{{\sin }^2}\theta \cos \theta d\theta } - 12\int_{\theta = 0}^{2\pi } {{{\cos }^2}\theta d\theta } + 36\int_{\theta = 0}^{2\pi } {{{\sin }^2}\theta d\theta + 2\int_{\theta = 0}^{2\pi } {\cos \theta d\theta } } } \right] \cr
& \,\,\, = 6\left[ { - 36\left( 0 \right) - 12\left( \pi \right) + 36\left( \pi \right) + 2\left( 0 \right)} \right] \cr
& \,\,\, = 144\pi \cr} $$
Next, compute the area element vector $d\bf{S}$ and $\nabla \times {\bf{F}}$
$$\eqalign{
& {\bf{x}} = x{\bf{i}} + y{\bf{j}} + \left( {1 - x} \right){\bf{k}} \cr
& d{\bf{S}} = \left( {{{\partial {\bf{x}}} \over {\partial x}} \times {{\partial {\bf{x}}} \over {\partial y}}} \right)dxdy = \left| {\matrix{
{\bf{i}} & {\bf{j}} & {\bf{k}} \cr
1 & 0 & { - 1} \cr
0 & 1 & 0 \cr
} } \right|dxdy = \left( {{\bf{i}} + {\bf{k}}} \right)dxdy \cr
& dS = \left\| {d{\bf{S}}} \right\| = \sqrt 2 dxdy \cr
& {\bf{n}} = {1 \over {\sqrt 2 }}\left( {{\bf{i}} + {\bf{k}}} \right) \cr
& \nabla \times {\bf{F}} = \left| {\matrix{
{\bf{i}} & {\bf{j}} & {\bf{k}} \cr
{{\partial _x}} & {{\partial _y}} & {{\partial _z}} \cr
{xy} & {2z} & {6y} \cr
} } \right| = 4{\bf{i}} - x{\bf{k}} \cr} $$
I think you had a mistake in this part $d{\bf{S}}=dS {\bf{n}}$ where $\sqrt2$ cancels. Finally, the surface integral will be
$$\eqalign{
& I = \int\!\!\!\int {\nabla \times {\bf{F}} \cdot d{\bf{S}}} = \int_{x = - 6}^6 {\int_{y = - \sqrt {36 - {x^2}} }^{\sqrt {36 - {x^2}} } {\left( {4 - x} \right)dydx} } \cr
& \,\,\,\, = \int_{x = - 6}^6 {2\left( {4 - x} \right)\sqrt {36 - {x^2}} dx} \cr
& \,\,\,\, = \int_{x = - 6}^6 {8\sqrt {36 - {x^2}} dx} = 8\int_{x = - 6}^6 {\sqrt {36 - {x^2}} dx} \cr
& \,\,\,\, = 8\left( {18\pi } \right) = 144\pi \cr} $$
Best Answer
Yes, you're right to want to use Stokes's Theorem. If $C$ is encloses the region $R$ in the plane $2x+3y+z=7$, with outward normal $\vec n = (2,3,1)/\sqrt{14}$, then $$\int_C \vec F\cdot d\vec r = \iint_R (\text{curl } \vec F\cdot\vec n)dA = \frac{18}{\sqrt{14}}\text{area}(R).$$ Can this ever have a maximum value?