# 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[55]:= ListPlot3D[z_list]
During evaluation of In[55]:= 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[55]:= ListPlot3D::arrayerr: z_list must be a valid array or a list of valid arrays. >>

Out[55]= 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". Dec 21, 2014 at 0:39
• You shouldn't use _ in parameter names. What happens if you replace your z_list with, e.g., zList? Dec 21, 2014 at 0:45
• Or to interactively vary the parameters and see the effects, use Manipulate. reference.wolfram.com/language/ref/Manipulate.html Dec 21, 2014 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. Dec 21, 2014 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? Dec 21, 2014 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[].