Plot3D list elements

Question

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]
    ]
   );

Answer 1

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, θ]}]
]