How Do I combine NDSolve (RungeKutta4thOrder) and Loop(For or Do) together?

thaikee84 · January 13, 2015, 7:15am

Hello everyone, previously I had ask about NDSolve with specified Runge Kutta 4th Order method for my equation. And here is what I got for combining few equations into one.

Z = 1

In[199]:= Clear[Z]

In[197]:= 
ClassicalRungeKuttaCoefficients[4, prec_] := 
 With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
   bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
  N[{amat, bvec, cvec}, prec]]


In[208]:= NDSolve[{Derivative[1][y][x] == 
   Piecewise[{{(y[x] + x^3 + 3 Z),
      0 <= x <= 1},
     
     {(y[x] + x^2 + 2 Z),
      1 <= x <= 2},
     
     {(y[x] + x + Z),
      2 <= x <= 3}}]
  , y[0] == 0}, y, {x, 0, 3}, 
 Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
   "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
 "StartingStepSize" -> 1/10]

Out[208]= {{y -> InterpolatingFunction[{{0., 3.}}, <>]}}

Z is a variable that I can change depending of the value for y[3] that I need.
Here’s the questions;

What command that I need to key-in in order for the equations to automatically generate the Plot and the value for y[3]?
As what I have been currently doing is I keep on clicking “get solution”, “plot” and " evaluate at point" buttons provided by Mathematica itself.
I need to use Loop command to find the value of Z so that I can get the value of y[3]=100 exactly. Currently I’m just try the value of Z randomly until I obtain y[3]=100 exactly. The For loop value of Z increment is 0.01 manually, the tutorials given to me was i++ command but the increment value is 1, how do I change it?

Thank you for your help guys.

thaikee84 · January 14, 2015, 6:08am

Hello,
I tried Break[] to get the value of z, when y[x] achieved 100,

For[z = 0,
  z < 100,
  z++,
  NDSolve[{Derivative[1][y][x] == y[x] + x^3 + 3 z, y[0] == 0}, 
   y, {x, 0, 3}]] If[y[3] > 100, Break[]]; y[3]

But the output is just y[3].

iamorozov1 · January 14, 2015, 12:54pm

Hi Thai,

Instead of loops try ParametricNDSolve[].

f = ParametricNDSolveValue[
	{
		Derivative[1][y][x] == 
		Piecewise[
			{
				{(y[x] + x^3 + 3 z),0 <= x <= 1},
				{(y[x] + x^2 + 2 z),1 <= x <= 2},
				{(y[x] + x + z),2 <= x <=3}
			}
		],
		y[0]==0
	},
	y[3.],
	{x,0.,3.},
	z
] ;

point = {z /. FindRoot[f[z] == 100.,{z,1},Evaluated->False], 100.} ;
Show[
	ListPlot[Evaluate[{#,f[#]} &/@ Range[1.,3.,0.05]],Joined->True,Frame->True,FrameLabel->{"Z","Y[3]"}],
	Graphics[{PointSize[Large],Red,Point[point]}],
	ImageSize -> 400
]

I.M.

thaikee84 · January 15, 2015, 1:41am

Thank you, Ivan, for your help.
This will definitely help me a lot in more than one way.

Is there a command that show the value of Z?, with the accuracy of my choice? for example I need to have the value of Z with the accuracy of 0.1 or 0.01 or 0.001, what command do I key-in?

By the way, I’m new to Mathematica, so I will never have think of this method.

Just to make confirmation, is ParametricNDSolveValue[] command, also use Runge-Kutta 4th Order Method as the base of the solution?

iamorozov1 · January 15, 2015, 6:18am

Just to make confirmation, is ParametricNDSolveValue[] command, also use Runge-Kutta 4th Order Method as the base of the solution?

No, but setting method is identical to that of NDSolve[], also do you need a fixed step rk4?

ClassicalRungeKuttaCoefficients[4,prec_] := With[
	{
		amat={{1/2},{0,1/2},{0,0,1}},
		bvec={1/6,1/3,1/3,1/6},
		cvec={1/2,1/2,1}
	},
		N[{amat,bvec,cvec},prec]
] ;

f = ParametricNDSolveValue[
	{
		Derivative[1][y][x] == 
		Piecewise[
			{
				{(y[x] + x^3 + 3 z),0 <= x <= 1},
				{(y[x] + x^2 + 2 z),1 <= x <= 2},
				{(y[x] + x + z),2 <= x <=3}
			}
		],
		y[0]==0
	},
	y[3.],
	{x,0.,3.},
	z,
	Method -> {
		"ExplicitRungeKutta",
		"DifferenceOrder" -> 4, 
		"Coefficients" -> ClassicalRungeKuttaCoefficients
	}, 
	StartingStepSize ->  1/10
];

Is there a command that show the value of Z?

The defined f functions takes value of Z and returns y[3] as output. Than the code:

FindRoot[f[z] == 100., {z, 1}, Evaluated -> False]

finds the value of Z for which y[3] = 100.

This creates a list {{z1,y3[z1]},{z2,y3[z2]},...} (see documentation for Range[] function)

{#, f[#]} & /@ Range[1., 3., 0.05]

Note, f can be plotted as a regular function of Z:

Plot[f[z], {z, 1, 3}]

I.M.

iamorozov1 · January 15, 2015, 6:33am

Also, you can modify f to return y instead of y[3.] (just remove [3.] in def. of f):

(* f[z][x] *)
f[1.5][3.]
f[#][3.] & /@ Range[1., 2., 0.1]
Plot[f[1.5][x], {x, 0, 3}]
ContourPlot[f[z][x], {x, 0, 3}, {z, 1, 3}]

I.M.

thaikee84 · January 15, 2015, 6:46am

Thanks again Ivan.
Your guidance really help me a lot.
I’m a student taking masters right now, so I need to confirm with the methods that I am using to solve the equations.
For my method of solution, I have determined to use RK4 with certain fixed step in my work, to achieve high accuracy.

iamorozov1 · January 15, 2015, 8:30am

To be sure that integration step is fixed you need to set “FixedStep” method as well as MaxStepFraction:

ClassicalRungeKuttaCoefficients[4,prec_] := With[
    {
       amat={{1/2},{0,1/2},{0,0,1}},
       bvec={1/6,1/3,1/3,1/6},
       cvec={1/2,1/2,1}
    },
       N[{amat,bvec,cvec},prec]
] ;

f = ParametricNDSolveValue[
    {
       Derivative[1][y][x] == 
       Piecewise[
         {
          {(y[x] + x^3 + 3 z), 0 <= x < 1},
          {(y[x] + x^2 + 2 z), 1 <= x < 2},
          {(y[x] + x + z), 2 <= x <= 3}
         }
       ],
       y[0] == 0
    },
    y,
    {x,0.,3.},
    z,
    Method -> {
		"FixedStep",
		Method -> {
			"ExplicitRungeKutta",
			"DifferenceOrder" -> 4, 
			"Coefficients" -> ClassicalRungeKuttaCoefficients
		}
	}, 
    StartingStepSize ->  #,
	MaxStepFraction -> #
] & ;

With this definition f is a function of step, z and x:

(* f[step][z][x] *)
f[0.5][1.][3.]
f[0.1][1.][3.]

You can then check that the step size is indeed fixed by requesting the grid in x used for integration:

f[0.5][1.]["Grid"]
f[0.1][1.]["Grid"]
(* commented output
{{0.}, {0.5}, {1.}, {1.}, {1.5}, {2.}, {2.}, {2.5}, {3.}, {3.}}    
{{0.}, {0.1}, {0.2}, {0.3}, {0.4}, {0.5}, {0.6}, {0.7}, {0.8}, {0.9}, {1.}, {1.}, {1.1}, {1.2}, {1.3}, {1.4}, {1.5}, {1.6}, {1.7}, {1.8}, {1.9}, {2.}, {2.}, {2.1}, {2.2}, {2.3}, {2.4}, {2.5}, {2.6}, {2.7}, {2.8}, {2.9}, {3.}, {3.}}  *)

Also note that due to Piecewise[] there are almost identical points (auto discontinuity handling) on x-grid, like {1.}, {1.} (in fact {0.9999999999999716},{1.}).

I.M.

thaikee84 · January 19, 2015, 1:43am

Hello Ivan,

I have encounter this error when I adjusted the starting point of my x0 value into higher or lower value.
Is it not, when I have increase the starting value, I would not encounter the Maximum Steps error?
I have tried few values and found the correct value that do not cause this error.
But, even though I had this error message, the value of Z root can still be found.

FindRoot[f,{x,Subscript[x, 0]}]

searches for a numerical root of f, starting from the point x=Subscript[x, 0].

FindRoot[f[z] == 100., {z, 1}, Evaluated → False]

ParametricNDSolveValue::mxst: Maximum number of 10000 steps reached at the point x$6275 == 0.6709750429824323`.

thaikee84 · January 30, 2015, 3:15am

Hello again,

If I have found the value of z

f = ParametricNDSolveValue[{Derivative[1][y][x] ==
Piecewise[{{(y + x^3 + 3 z),
0 <= x <= 1}, {(y + x^2 + 2 z),
1 <= x <= 2}, {(y + x + z), 2 <= x <= 3}}], y[0] == 0},
y[3.], {x, 0., 3.}, z];
FindRoot[f[z] == 100., {z, 1}, Evaluated → False]

And Mathematica shows the answer in this form

{z → 1.7049}

How do I use the value of the z to continue my equation?
For example, the equation is
A=100-z

thaikee84 · January 30, 2015, 7:02am

Hello again,

For this equation, I would like to do an iteration to increase Twater temperature until the Tnew are the same.
I would like to know which type of loop that I can use until Twater=Tnew
Shall I use ParametricNDSolve or other type of method?

(*Temperature of water*)
Twater = 300

ClassicalRungeKuttaCoefficients[4, prec_] := With[
   {
    amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
    bvec = {1/6, 1/3, 1/3, 1/6},
    cvec = {1/2, 1/2, 1}
    },
   N[{amat, bvec, cvec}, prec]
   ];
f = ParametricNDSolveValue[{Derivative[1][y][x] ==
     Piecewise[{
       {(y[x] + x^3 + 3 z - Twater), 0 <= x <= 1},
       {(y[x] + x^2 + 2 z), 1 <= x <= 2},
       {(y[x] + x + z), 2 <= x <= 3}
       }
      ]
    , y[0] == 0
    },
   y[3.],
   {x, 0., 3.},
   z,
   Method -> {"ExplicitRungeKutta",
     "DifferenceOrder" -> 4, 
     "Coefficients" -> ClassicalRungeKuttaCoefficients
     },
   StartingStepSize -> 1/10
   ];
point = {z /. FindRoot[f[z] == 100., {z, 1}, Evaluated -> False], 
   100.};
Show[ListPlot[Evaluate[{#, f[#]} & /@ Range[60., 100., 0.05]], 
  Joined -> True, Frame -> True, FrameLabel -> {"Z", "Y[3]"}], 
 Graphics[{PointSize[Large], Red, Point[point]}], ImageSize -> 400]

FindRoot[f[z] == 100., {z, 1}, Evaluated -> False]

Tnew = 170 + 2*z /. FindRoot[f[z] == 100., {z, 1}, Evaluated -> False]

iamorozov1 · January 30, 2015, 7:36am

Hi,

you can add Twater as a parameter into def. of f :

(*Temperature of water*)
Twater = 300
(*Dimensions of the MHP*)
ClassicalRungeKuttaCoefficients[4, prec_] := With[
   {
    amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
    bvec = {1/6, 1/3, 1/3, 1/6},
    cvec = {1/2, 1/2, 1}
    },
   N[{amat, bvec, cvec}, prec]
   ];
f = ParametricNDSolveValue[{Derivative[1][y][x] ==
     Piecewise[{
       {(y[x] + x^3 + 3 z - t), 0 <= x <= 1},
       {(y[x] + x^2 + 2 z), 1 <= x <= 2},
       {(y[x] + x + z), 2 <= x <= 3}
       }
      ]
    , y[0] == 0
    },
   y[3.],
   {x, 0., 3.},
   {z,t},
   Method -> {"ExplicitRungeKutta",
     "DifferenceOrder" -> 4, 
     "Coefficients" -> ClassicalRungeKuttaCoefficients
     },
   StartingStepSize -> 1/10
   ];
(* Tnew *)
170 + 2*z /. FindRoot[f[z,Twater] == 100., {z, 1}, Evaluated -> False]
(* Find eq. temp. *)
g = 170 + 2 z /. FindRoot[f[z,#] == 100., {z, 1}, Evaluated -> False] & ;
Teq = x /. FindRoot[g[x]==x,{x,Twater},Evaluated -> False]
(* Check *)
170 + 2*z /. FindRoot[f[z,Teq] == 100., {z, 1}, Evaluated -> False]

I.M.

thaikee84 · January 30, 2015, 7:58am

Thanks again, Ivan.

May I know how does the equation work?
At first, the equation help me solve z with initial Twater = 300

f = ParametricNDSolveValue[{Derivative[1][y][x] ==
     Piecewise[{
       {(y[x] + x^3 + 3 z - t), 0 <= x <= 1},
       {(y[x] + x^2 + 2 z), 1 <= x <= 2},
       {(y[x] + x + z), 2 <= x <= 3}
       }
      ]
    , y[0] == 0
    },
   y[3.],
   {x, 0., 3.},
   {z,t},
   Method -> {"ExplicitRungeKutta",
     "DifferenceOrder" -> 4, 
     "Coefficients" -> ClassicalRungeKuttaCoefficients
     },
   StartingStepSize -> 1/10
   ];

Then it help me to look for a Tnew

(*Tnew*)
170 + 2*z /. 
 FindRoot[f[z, Twater] == 100., {z, 1}, Evaluated -> False]

With the value is not True,
The third function repeat ParametricNDSolveValue until it reaches True

(* Find eq. temp. *)
g = 170 + 2 z /. FindRoot[f[z,#] == 100., {z, 1}, Evaluated -> False] & ;
Teq = x /. FindRoot[g[x]==x,{x,Twater},Evaluated -> False]

Then finally, as it reaches True,
Recheck the Twater and conclude the correct value for Twater= Tnew is True

(* Check )
170 + 2z /. FindRoot[f[z,Teq] == 100., {z, 1}, Evaluated → False]

Is that how it works?
Btw, what does # represents?

thaikee84 · January 30, 2015, 8:21am

So, here’s a little changes.
For my equation, there’s some information that I need to extract from Mathematica accordingly to the temperature of the water.
For example, as the water temperature changes, the thermal conductivity of water will change as well .

I received error because Twater is not a value yet and I’m stuck.

Then again, sorry for the previous question that I did not explain fully for what I was trying to do.

 (*Temperature of water*)
     Twater = 300
     Cwater = QuantityMagnitude[
       ThermodynamicData["Water", 
        "ThermalConductivity", {"Temperature" -> 
          Quantity[Twater, "Kelvins"]}]]
 
 ClassicalRungeKuttaCoefficients[4, prec_] := 
   With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
     bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
    N[{amat, bvec, cvec}, prec]];
 f = ParametricNDSolveValue[{Derivative[1][y][x] == 
      Piecewise[{{(y[x] + Cwater/3 x^3 + 3 z ), 
         0 <= x <= 1}, {(y[x] + Cwater/3 x^2 + 2 z), 
         1 <= x <= 2}, {(y[x] + Cwater/3 x + z), 2 <= x <= 3}}], 
     y[0] == 0}, y[3.], {x, 0., 3.}, {z, Cwater}, 
    Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
      "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
    StartingStepSize -> 1/10];
 
 FindRoot[f[z, Twater] == 100., {z, 1}, Evaluated -> False]
 (*Tnew*)
 170 + 2*z /. FindRoot[f[z, Twater] == 100., {z, 1}, Evaluated -> False]
 (*Find eq.temp.*)
 g = 170 + 2 z /. 
     FindRoot[f[z, #] == 100., {z, 1}, Evaluated -> False] &;
 Teq = x /. FindRoot[g[x] == x, {x, Twater}, Evaluated -> False]
 (*Check*)
 170 + 2*z /. FindRoot[f[z, Teq] == 100., {z, 1}, Evaluated -> False]

iamorozov1 · January 30, 2015, 8:31am

Btw, what does # represents?

# is a slot, it is used to define pure functions, e.g. :

f1[x_] = x^2 ;
f2 = #^2 & ;
(* the returned result is the same *)
{f1[2],f2[2]}
(* internal representation is different *)
DownValues[f1]
DownValues[f2]

See this post for more info as well as documentation.

Here is a cleaner (I hope) version of your code:

ClassicalRungeKuttaCoefficients[4, prec_] := With[
   {
    amat = {{1/2}, {0, 1/2}, {0, 0, 1}},
    bvec = {1/6, 1/3, 1/3, 1/6},
    cvec = {1/2, 1/2, 1}
    },
   N[{amat, bvec, cvec}, prec]
] ;

(* Define function f : {z,t} ->  y[3.] *)
f = ParametricNDSolveValue[
	{
	Derivative[1][y][x] ==
     Piecewise[{
       {(y[x] + x^3 + 3 z - t), 0 <= x <= 1},
       {(y[x] + x^2 + 2 z), 1 <= x <= 2},
       {(y[x] + x + z), 2 <= x <= 3}
       }
      ]
    , y[0] == 0
    },
   y[3.],
   {x, 0., 3.},
   {z,t},
   Method -> {"ExplicitRungeKutta",
     "DifferenceOrder" -> 4, 
     "Coefficients" -> ClassicalRungeKuttaCoefficients
     },
   StartingStepSize -> 1/10
   ];

(* Define function g : t_initial -> t_final *)
g = 170 + 2 z /. FindRoot[f[z,#] == 100., {z, 1}, Evaluated -> False] & ;

(* Find fix point of g, i.e. g[T] = T *)
Tinitial = 300. ;
Tfinal = x /. FindRoot[g[x] == x,{x,Tinitial},Evaluated -> False]
(* Or *)
Quiet[FixedPoint[g,Tinitial]]

I.M.

thaikee84 · January 30, 2015, 8:42am

Thank you very much.

iamorozov1 · January 30, 2015, 9:10am

You can define Cwater as a function of temp. :

    Cwater = QuantityMagnitude[
    	ThermodynamicData[
    		"Water",
    		"ThermalConductivity", {"Temperature" ->Quantity[#, "Kelvins"]}
    	]
    ] & ;
    (* eg: *)
    Cwater[300.]
    Cwater[310.]
    Cwater[320.]

Then f : {z,c} -> y[3] where c is a parameter (former Cwater) can be called as :

f[1.9,Cwater[300.]]

I.M.

thaikee84 · February 2, 2015, 3:22am

Hello,

Please forgive me because I’ve tried some of the methods but it still doesn’t work for me. Or, I just misunderstand some of the explanations

This is the original code that was written and I receive the error for Cwater = 0.61… cannot be used as a parameter because Cwater gain a fixed value.

   (*Temperature of water*)
         Twater = 300
         Cwater = QuantityMagnitude[
           ThermodynamicData["Water", 
            "ThermalConductivity", {"Temperature" -> 
              Quantity[Twater, "Kelvins"]}]]
    
 ClassicalRungeKuttaCoefficients[4, prec_] := 
   With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
     bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
    N[{amat, bvec, cvec}, prec]];
 f = ParametricNDSolveValue[{Derivative[1][y][x] == 
      Piecewise[{{(y[x] + Cwater/3 x^3 + 3 z ), 
         0 <= x <= 1}, {(y[x] + Cwater/3 x^2 + 2 z), 
         1 <= x <= 2}, {(y[x] + Cwater/3 x + z), 2 <= x <= 3}}], 
     y[0] == 0}, y[3.], {x, 0., 3.}, {z, Cwater}, 
    Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
      "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
    StartingStepSize -> 1/10];

 FindRoot[f[z, Twater] == 100., {z, 1}, Evaluated -> False]
 (*Tnew*)
 170 + 2*z /. FindRoot[f[z, Twater] == 100., {z, 1}, Evaluated -> False]
 (*Find eq.temp.*)
 g = 170 + 2 z /. 
     FindRoot[f[z, #] == 100., {z, 1}, Evaluated -> False] &;
 Teq = x /. FindRoot[g[x] == x, {x, Twater}, Evaluated -> False]
 (*Check*)
 170 + 2*z /. FindRoot[f[z, Teq] == 100., {z, 1}, Evaluated -> False]

Then I changed the command for Cwater by

Cwater = QuantityMagnitude[
    ThermodynamicData["Water", 
     "ThermalConductivity", {"Temperature" -> 
       Quantity[#, "Kelvins"]}]] &;
(*eg:*)
Temp = Quantity[{t, 300, 360}, "Kelvins"]

And I received this error

ThermodynamicData::quants: List {t,300,360} does not consist of real numbers. >>

Shall I change into For loop with Temperature as i++ until 400 and stop when Tinitial = Tfinal?
How do I do that?

thaikee84 · February 3, 2015, 9:06am

Hello again,

Now I have tried this method manually, it’s suppose to work. Just wondering how do I do it automatically. For example

(*Temperature of water*)
Twater = 500
Cwater = QuantityMagnitude[
  ThermodynamicData["Water", 
   "ThermalConductivity", {"Temperature" -> 
     Quantity[Twater, "Kelvins"]}]]

ClassicalRungeKuttaCoefficients[4, prec_] := 
  With[{amat = {{1/2}, {0, 1/2}, {0, 0, 1}}, 
    bvec = {1/6, 1/3, 1/3, 1/6}, cvec = {1/2, 1/2, 1}}, 
   N[{amat, bvec, cvec}, prec]];

f = ParametricNDSolveValue[{Derivative[1][y][x] == 
     Piecewise[{{(y[x] + x^3 + 3 z - 18000*Cwater), 
        0 <= x <= 1}, {(y[x] + x^2 + 2 z), 
        1 <= x <= 2}, {(y[x] + x + z), 2 <= x <= 3}}], y[0] == 0}, 
   y[3.], {x, 0., 3.}, z, 
   Method -> {"ExplicitRungeKutta", "DifferenceOrder" -> 4, 
     "Coefficients" -> ClassicalRungeKuttaCoefficients}, 
   StartingStepSize -> 1/10];
point = {z /. FindRoot[f[z] == 100., {z, 1}, Evaluated -> False], 
   100.};
Show[ListPlot[Evaluate[{#, f[#]} & /@ Range[60., 100., 0.05]], 
  Joined -> True, Frame -> True, FrameLabel -> {"Z", "Y[3]"}], 
 Graphics[{PointSize[Large], Red, Point[point]}], ImageSize -> 400]

FindRoot[f[z] == 100., {z, 1}, Evaluated -> False]

Tnew = 170 + 2*z /. FindRoot[f[z] == 100., {z, 1}, Evaluated -> False]

How do I do this command?:

Tinitial.
Found Tnew.
Tinitial = Tnew, False.
Replace Tinitial with Tnew
Repeat.
Tinitial = Tnew, True.
Stop.
OR Stop when iteration is 100 and no suitable value is found

NOTE: Tried this equation with few iterations, Tnew keep getting further away from Tinitial.

brucem · February 3, 2015, 6:22pm

Similar discussion in http://community.wolfram.com/groups/-/m/t/432065 .

Topic		Replies	Views
How do I repeat/loop multiple functions that are related? Questions engineering , mathematica , wolfram-alpha , calculus , equation-solving , optimization , numerical-computation	12	34970	March 9, 2015
Runge Kutta method to solve the Equation Questions calculus , equation-solving , numerical-computation	7	15540	February 5, 2015
ParametricNDSolveValue error Questions ruliology , computer-science , engineering , physics , mathematica , calculus , numerical-computation	0	9772	March 2, 2015
Solve a system of differential equations in a loop Questions physics , calculus , equation-solving , numerical-computation	7	5592	March 23, 2021
Caputo-Hadamard fractional derivatine nonlinear equation system Questions calculus , equation-solving , graphics-and-visualization , numerical-computation	1	2583	January 10, 2023

How Do I combine NDSolve (RungeKutta4thOrder) and Loop(For or Do) together?

Related topics