# How to solve correctly inside a DO loop and then plot

Here is the corresponding code

First the equations setup

Vcl[x_, y_, z_] := (-G*(Mcl/a))/Sqrt[(x^2 + y^2 + z^2)/a^2 + c^2];
V[x_, y_, z_] := Vcl[x, y, z] + 1/2*(κ2 - 4*ω^2)*x^2 + 1/2*v2*z^2;

G = 1;
Mcl = 2.2; a = 0.182; c = 1;
ω = 1; κ2 = 1.8; v2 = 7.6;
cl = -3.264444506;
ch = 0.1;
E0 = cl*(1 - ch);


Then the DO loop with the Solve

data = {};
Do[
Vx = V[x, 0, z0];
sol = Solve[Vx == E0 && -1.01 < x < 1.01,x];
xmin = x /. sol[[1]];
xmax = x /. sol[[2]];
AppendTo[data, {xmin, xmax, z0}],
{z0, 0, 0.6, 0.001}
]


And finally the plot

neg = data[[All, {1, 3}]];
pos = data[[All, {2, 3}]];
L1 = ListPlot[neg, Joined -> True,
PlotStyle -> {Black, Thickness[0.003]}];
L2 = ListPlot[pos, Joined -> True,
PlotStyle -> {Black, Thickness[0.003]}];
L0 = Show[{L1, L2}, Axes -> False, Frame -> True,
PlotRange -> {{-1, 1}, {0, 0.6}}, ImageSize -> 550]


which produces this

The plot has several issues: (i) it should be a continuous line from -1 to 1, while we observe a gap near 0 and also it never reaches +1, (ii) the lower horizontal line it should not be there; probably Joined joins two extreme points. I suspect, that both issues are related to the Solve inside the DO loop. There z0 belongs to the interval [0, 0.6] but the equation does not have solutions for all these values. The data list should somehow store only the correct solutions and reject the cases where there are no solutions at all.

• There are some warnings when solving your equation in the Do loop. You'd better correct the code. By the way, the solutions to your equations contains complex numbers, so the condition -1.01 < x < 1.01 might give you the void set. Apr 10 '14 at 8:28
• @Kuba I think that the flood of error messages is because for some values of z0 there are no solutions. This is exactly what I want to know; how to get reject these cases from the list. The Solve syntax is corrected. BTW, thanks for the shorter version of neg! Apr 10 '14 at 8:36
• @Z-Y.L I think that the warnings when solving the equations in the Do loop are mainly because there are no solutions for every value of z0; Apr 10 '14 at 8:38
• @Kuba So is there no solution to the problem at all? Apr 10 '14 at 8:54
• Why not use ContourPlot[V[x, 0, z0] == E0, {x, -2, 2}, {z0, 0, 0.6}] Apr 10 '14 at 10:05

You could avoid the problem altogether by using ContourPlot:

ContourPlot[V[x, 0, z0] == E0, {x, -2, 2}, {z0, 0, 0.6}]


Using a finer grid in the loop (e.g. 0.0001) helps to reduce the gap in the middle.

If we do not join the data points and give different colors to your neg and pos lists, we see that something strange happens for small $z_0$: Your equation has 4 solutions (two negative, two positive), so just using the first or the second argument wont pick the smallest or largest.

 neg = data[[All, {1, 3}]];
pos = data[[All, {2, 3}]];
L1 = ListPlot[neg, PlotStyle -> {Blue, Thickness[0.002]}];
L2 = ListPlot[pos, PlotStyle -> {Red, Thickness[0.002]}];
L0 = Show[{L1, L2}, Axes -> False, Frame -> True,
PlotRange -> {{-1, 1}, {0, 0.6}}, ImageSize -> 550]


We can also see this if we plot $V$ and $E_0$ together for different $z_0$ (just at the beginning you can see it)

 Export["~/test.gif", Table[Plot[{V[x, 0, z0], E0}, {x, -1.01, 1.01}, PlotRange -> {-4.5, -1.5}, Frame -> True], {z0, 0.2, 0.55, 0.01}], "GIF"]


If you want xmin and xmax you will always have a gap close to +/- 1, but you can do it like this:

 data = {};
Do[Vx = V[x, 0, z0];
sol = Solve[Vx == E0 && -1.01 < x < 1.01, x];
xmin = x /. First@sol;
xmax = x /. Last@sol;
AppendTo[data, {xmin, xmax, z0}], {z0, 0.25, 0.55, 0.0001}]


If you want all real solutions between -1.01 and 1.01 you can do this:

data = {};
Do[Vx = V[x, 0, z0];
sol = Solve[Vx == E0 && -1.01 < x < 1.01, x];
allx = {x, z0} /. sol;
If[NumberQ[First@Flatten@allx],
AppendTo[data, Flatten[{allx}, 1]]],
{z0, 0.25, 0.55, 0.0001}]
ListPlot[Flatten[data, 1]]


where it seems a little hacky with all the Flatten and I am checking if there is a solution before I append it to the list, but it works and you get (maybe) what you want (where using a finer grid for z0 reduces the gap even more):

• This would be improved if you use Reap/Sow instead of AppendTo Apr 10 '14 at 18:21
• sure, I just didn't want to obfuscate the important changes, that help to solve the main problem. Apr 11 '14 at 5:22

Here is a 'cheat':

p1 = Plot3D[V[x, 0, z], {x, -1.01, 1.01}, {z, 0, 0.6},
MeshFunctions -> (#3 &), Mesh -> {{E0}}, MeshStyle -> {Red, Thick},
PlotPoints -> 100]


You can then extract the desired mesh points from the object:

lns = Cases[p1, Line[x_] :> x, Infinity];
gr = Graphics[{Red, Thick,
Line[(p1[[1, 1]][[lns[[2]]]])[[All, {1, 2}]]]}, Axes -> True,
Frame -> True, AspectRatio -> 1, AxesOrigin -> {0, 0}]


• See above a much simpler solution by @Simon Woods! Apr 10 '14 at 10:25
• @Vaggelis_Z thank you...should have thought of that...and should have read the comments... Apr 10 '14 at 10:28
• @ubpdqn at least you're keeping it funky :) Apr 10 '14 at 16:18
• Ingenuous yet horribly convoluted. I like it! +1 Apr 10 '14 at 17:12