# Plot3D list elements

I've been having some trouble with Plot3D. I have a quite complex function of two variables (actually three but the first one doesn't matter) that returns a list with two elements. Then I wrote the following code to see the behaviour of the first element:

Module[{R = 1},
Plot3D[
g[R, θ, z][[1]],
{θ, 0, Pi}, {z, 0, 4},
AxesLabel -> {θ, z, Subscript[F, θ]}
]
]


Well the thing is: when the coordinate $\theta=0$, the first element is zero and the plot doesn't show that:

I know that the problem isn't from the function because if I plot the first element keeping $\theta=0$:

Module[{R = 1},
Plot3D[
g[R, 0, z][[1]],
{θ, 0, Pi}, {z, 0, 4},
AxesLabel -> {θ, z, Subscript[F, θ]}
]
]


(in the code I only changed that $\theta \to 0$) it outputs correctly:

Moreover, If I use Manipulate it behaves as expected; so I must be doing something wrong with Plot3D... can anyone help me?

PS: Here is the horrible definition of the function $g$:

g[R_, \[Theta]_,
z_] := (Module[{For1 = 0, For2 = 0, d1 = 0, d2 = 0, Forca = {0, 0},
nmax = 10},

If[\[Theta] == 0,
For[n = 0, n < nmax,
d1 = Sqrt[R^2*(2 Pi*n)^2 + z^2];
If[n == 0,
If[d1 + 0.0 < R*Pi,
For1 = -1/(2 Pi*d1^2);
Vec1 = {0, z};
,
For1 = -1/(4*d1^2*ArcSin[R*Pi/d1]);
Vec1 = {0, z};
];
,
For1 = -1/(4*d1^2*ArcSin[R*Pi/d1]);
Vec1 = {0, z};
For2 = -1/(4*d1^2*ArcSin[R*Pi/d1]);
Vec2 = {0, z};
];
Forca += For1*Vec1 + For2*Vec2;
++n;
];
(*Here works iff \[Theta]=0*)
,
For[n = 0, n < nmax,
d1 = Sqrt[R^2*(\[Theta] + 2 Pi*n)^2 + z^2];
d2 = Sqrt[R^2*(-2 Pi*(n + 1) + \[Theta])^2 + z^2];

If[d1 + 0.0 < R*Pi,
For1 = -1/(2 Pi*d1^2);
Vec1 = {\[Theta], z};
,
For1 = -1/(4*d1^2*ArcSin[R*Pi/d1]);
Vec1 = {\[Theta] + 2 Pi*n, z};
];

If[d2 + 0.0 < R*Pi,
For2 = -1/(2 Pi*d2^2);
Vec2 = {-2 Pi*(n + 1) + \[Theta], z};
,
For2 = -1/(4*d2^2*ArcSin[R*Pi/d2]);
Vec2 = {-2 Pi*(n + 1) + \[Theta], z};
];
Forca += For1*Vec1 + For2*Vec2;
++n;
]
(*Here works for all points except \[Theta]=0*)
];
N[Forca]
]
);

-
@Jens I was answering a question in another comment, now deleted. Deleting –  belisarius Mar 10 '13 at 19:40

This is probably due to the degenerate nature of your function.

Take for example this slice (z = 1):

Plot[g[1, θ, 1][[1]], {θ, -2, 2}, PlotRange -> All]


It looks like it's continuous a θ=0 but it isn't:

g[1, #, 0.5][[1]] & /@ {-0.0001, 0, -0.0001}


{0.07960574141, 0., 0.07960574141}

Things gets worse for lower z values:

Plot[g[1, θ, 0.01][[1]], {θ, -2, 2}, PlotRange -> All]


The zero value at θ=0 is probably hard to find with most samplings of the parameter space.

A workaround would be to use ListPlot3D:

Module[{R = 1},
ListPlot3D[
Flatten[Table[{θ, z, g[R, θ, z][[1]]}, {z, 0.001, 4, 0.1}, {θ, 0, Pi, Pi/40}], 1],
AxesLabel -> {θ, z, Subscript[F, θ]}]
]


-