# How to plot output data of Do-command?

The following in a simple numerical solution of the simple harmonic oscillator.

  \[Epsilon] = 0.05; x = 0; f = 0; g = 1; i = 0;
Do[
f1 = f + \[Epsilon]*g;
g1 = g - \[Epsilon]*f;
x1 = x + \[Epsilon]; i1 = i + 1;
Print[x1, " ", f1, " ", g1];
x = x1; f = f1; g = g1; i = i1,
{i, 0, 100}]


What I want to do is to plot the outputs. That is to have a graph with f1 on the y-axis and x on the x-axis and another one with g1 on the y-axis and x on the x-axis. I am not sure how to do this. I read in another post here someone suggesting to make a list of all the data and then plot them, but I am not sure I understand how to do this.

• None of the data {f1, g1, x}-values are saved anywhere. Once you finish the Do loop all you're left with are the last two values for each. Are you interested in creating a numerical solver from scratch (i.e,, is that why you're not using NDSolve)? Are you tied to using a Do loop for this? Or are you interested in the solution itself? – aardvark2012 Jun 24 '17 at 11:35
• I tried to use NDSolve, but the solution does not match the theory and the result I should get and I I am trying different ways of getting there. – DiSp0sablE_H3r0 Jun 26 '17 at 9:36

This is a good candidate for NestList

e = 0.05;
f = 0;
g = 1;
n = 100;

ListLinePlot[
Transpose@Most@NestList[{First@# + e Last@#, Last@# - e First@#} &, {f, g}, n],
DataRange -> {e, e*n},
PlotLegends -> {"f", "g"}]


• Thanks for your time. Your reply is as good as the one above. – DiSp0sablE_H3r0 Jun 24 '17 at 17:20
• my problem is the following. I have this [Epsilon] = 0.001; x = 0; f = 0.4; g = 0; i = 0; Do[ f1 = f + [Epsilon]*g; g1 = g + [Epsilon]*(-(2/(x + 0.1))*g + k[f]); x1 = x + [Epsilon]; i1 = i + 1; Print[x1, " ", f1, " ", g1]; x = x1; f = f1; g = g1; i = i1, {i, 0, 100000}] And then I used the command you wrote here to plot it. ListLinePlot[ Transpose@ Most@NestList[{First@# + e Last@#, Last@# - e First@#} &, {f, g}, n], DataRange -> {e, e*n}, PlotLegends -> {"f", "g"}] The problem. The data seem to converge to a solution, but the plot is oscillatory – DiSp0sablE_H3r0 Jun 26 '17 at 9:32
• What is k[f] ? I would pose this as a new question showing all steps – eldo Jun 26 '17 at 10:41
• k[f] is a very complicated expression. I will follow your suggestion and pose it as a new question. Thanks for the hint. – DiSp0sablE_H3r0 Jun 26 '17 at 11:05
• here is the link with the new question mathematica.stackexchange.com/questions/149109/… thank you once again for your time – DiSp0sablE_H3r0 Jun 26 '17 at 11:41

You could use AppendTo to collect the output and then plot. However, I think you have made your recursion unnecessarily complex. You could do something like this:

v = {0, 0, 1}
func[eps_] := {#1 + eps,
{#2, #3}.{1, eps},
{#3, #2}.{1, -eps}} &;


Then achieve your aim with NestList, e.g.:

Manipulate[dat = NestList[func[eps] @@ # &, v, 100];
Row[{TableForm[dat[[1 ;; 10]],
TableHeadings -> {None, {"x", "f", "g"}}],
ListPlot[{dat[[All, {1, 2}]], dat[[All, {1, 3}]]}, Joined -> True,
PlotLegends -> {"f", "g"}, ImageSize -> 200]}], {eps, {0.01, 0.05,
0.1, 0.2}}]


• Thank you for that. Very helpful and seems wonderful – DiSp0sablE_H3r0 Jun 24 '17 at 17:17

Your system of recursions can be solved exactly using RSolve

Clear[ϵ, f, g, x];

{f[ϵ_, n_], g[ϵ_, n_], x[ϵ_, n_]} =
{f[ϵ, n], g[ϵ, n], x[ϵ, n]} /. RSolve[{
f[ϵ, n] ==
f[ϵ, n - 1] + ϵ*g[ϵ, n - 1],
g[ϵ, n] ==
g[ϵ, n - 1] - ϵ*f[ϵ, n - 1],
x[ϵ, n] == x[ϵ, n - 1] + ϵ,
x[ϵ, 0] == 0, f[ϵ, 0] == 0, g[ϵ, 0] == 1},
{f[ϵ, n], g[ϵ, n], x[ϵ, n]}, n][[1]] //
ComplexExpand[#, TargetFunctions -> {Re, Im}] & //
FullSimplify

(*  {(1 + ϵ^2)^(n/2)*Sin[n*ArcTan[ϵ]],
(1 + ϵ^2)^(n/2)*Cos[n*ArcTan[ϵ]], n*ϵ}  *)


These expression can be evaluated for all real {n, ϵ}

Manipulate[
Row[{
Join[{Style[#, Bold, 12] & /@ {"n", "x", "f", "g"}},
Table[{n, x[ϵ, n], f[ϵ, n], g[ϵ, n]}, {n, 0,
9}]] // Grid[#, Alignment -> Left, Frame -> All] &,
Spacer[10],
Plot[{f[ϵ, n], g[ϵ, n]}, {n, 0, 130},
Frame -> True, Axes -> False,
FrameLabel -> {Style["n", Bold, 12], None},
PlotLegends -> {"f[ϵ,n]", "g[ϵ,n]"},
ImageSize -> 270]}],
{{ϵ, 0.05}, {0.01, 0.02, 0.05, 0.1, 0.2}},
ContentSize -> {600, 250}]


• Thank you for your time and your reply. I think it seems a bit more complicated than the two solutions above, but the plots seem great. – DiSp0sablE_H3r0 Jun 24 '17 at 17:21