# How do I plot results of complicated calculations in Mathematica?

I am new to Mathematica and am struggling with making plots of how features of solutions vary with parameters. For example, I can find the eigenvalues and zeroes of a Sturm Liouville problem (showing just the first nonzero eigenvalue for this example) with

sol = NDSolve[
{g''[x] -
4*c*(-((3 r sigma^3)/(L - x)^4) + (3 sigma^3)/(L + x)^4)*
g'[x] + mu[x]*g[x] == 0, mu'[x] == 0,
g[-b] == 1, g'[-b] == 0, g'[b] == 0},
{g, mu}, x,
Method -> {"EventLocator", "Event" -> g[x],
"EventAction" :>
Print["mu= ", lambda = mu[x], ",  g[", zero = x, "] = ",
g[x]]}
]


but I would like to see how these features vary with the parameters. Looking at the basic help files that come with Mathematica, I can loop over parameters using Do[(block of code given above),{c,cstart,cend,dc},{r,rstart,rend,dr}] etc. But then the documentation tells me how to make graphs by either typing thousands of numbers in curly braces or reading them from a file (of not obviously documented format).

I've spent two hours trying to store these numbers in arrays that ListPlot3d will recognize but I either get error messages when trying to create the arrays or error messages when ListPlot3d tries to read them. The closest I have come to success is writing the results to a sequential file

sfile = OpenWrite["zero.dat"]; Do[....
Write[sfile, {r, c, zero}]
....
Close[sfile]
z_list = Import["zero.dat", "Table"]


This gets the numbers in and out. If the imported data is displayed in table form it looks like a table of x,y,z data for ListPlot3d. But when I try to run it I get

In:= ListPlot3D[z_list]

During evaluation of In:= Rule::rhs: Pattern z_list appears on the right-hand side of rule SystemProtoPlotDumpmodelData$26641[WrappedValues,z_list,ChartingPrivateTag$26644]->z_list. >>

During evaluation of In:= ListPlot3D::arrayerr: z_list must be a valid array or a list of valid arrays. >>

Out= ListPlot3D[z_list]


The idea that I have to write results to a file to plot them sounds like I should be invoking WATFOR. There must be an elegant way to plot the results of complex calculations, but in two hours looking through help (and google), I haven't found it. How do I do this?

• What version of Mathematica do you have? Since version 9, there is ParametricNDSolve and WhenEvent has replaced "EventLocator". – Michael E2 Dec 21 '14 at 0:39
• You shouldn't use _ in parameter names. What happens if you replace your z_list with, e.g., zList? – Karsten 7. Dec 21 '14 at 0:45
• Or to interactively vary the parameters and see the effects, use Manipulate. reference.wolfram.com/language/ref/Manipulate.html – Bob Hanlon Dec 21 '14 at 0:54
• Your code is confusing, but at the very least you can simplify it by replacing the constraint mu'[x] == 0 with mu == mu0 (a constant). Don't compound problems by worrying about how to read and write variables (Import and Export). Give us the simplest version of your code that works. – David G. Stork Dec 21 '14 at 1:30
• Have you looked at the help file for NDSolve? There are many examples there that plot the output of the ODE. Maybe you could mimic one of these? – bill s Dec 21 '14 at 3:42

I woudl suggest defining your soulutions as a delayed function using SetDelayed (":=")

sol[c_, sigma_, L_, r_, b_] := ...


Now you can go ahead and generate a table of solutions like that:

sols = ParallelTable[sol[c, sigma, L, r, b], {c, cmin, cmax, dc}, {sigma, sigmamin, sigmamax, dsigma}, {L, Lmin, Lmax, dL}, {r, rmin, rmax, dr}, {b, bmin, bmax, db}]


and access a solution within sols by index such that

sols[[1,2,3,4,5]]


yields

sol[cmin + 0*dc, sigmamin + 1*dsigma, Lmin + 2*dL, rmin + 3*dr, bmin + 4*db]


you can also take subsets of the table you created and listplot them. I hope I understood your problem right.

PS: since your problem comes with 5 Parameters you get a 5 dimensional result wich is hard to plot, plotting subsets might be a solution. A good way to plot up to 3 dimensional results is using Image3D[].