# Get value from InterpolatingFunction

I'm trying to find My[t], Mx[t] and Mz[t] values when evaluated at t = Infinity.

Code to date:

γ = 2.675*10^8;
T1 = 100*10^-3;
T2 = 10*10^-3;
B0 = 2;
M0 = 1*10^-1;
B1 = 1*10^-7;

(*Above are constants*)

NDSolve[{Mx'[t] == -Mx[t]/T2, Mx[0] == 0}, Mx, {t, 0, 0.5}]

NDSolve[{My'[t] == γ*Mz[t]*B1 - My[t]/T2,
Mz'[t] == -γ*My[t]*B1 + (M0 - Mz[t])/T1,
My[0] == 0, Mz[0] == M0}, {My, Mz}, {t, 0, 0.5},
MaxSteps -> Infinity]

Plot[Evaluate[Mx[t] /. %], {t, 0, 0.5},
FrameLabel -> {Style[t, 14], Style[Subscript[M, x], 12]},
PlotTheme -> "Scientific"]

Plot[Evaluate[{Mz[t], My[t]} /. %], {t, 0, 0.5},
FrameLabel -> {Style[t, 14], Style[M, 14]},
PlotTheme -> "Scientific",
PlotLegends ->
Placed[{"\!$$\*SubscriptBox[\(M$$, $$y$$]\)",
"\!$$\*SubscriptBox[\(M$$, $$z$$]\)"}, {0.9, 0.85}],
ClippingStyle -> None, PlotRange -> All,
FrameTicks -> {{Automatic, All}, {Automatic, None}}]

• Welcome to Mathematica.SE! 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2) Take the tour and check the faqs! 3) When you see good questions and answers, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. Also, please remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign!
– user9660
Nov 6, 2015 at 16:20
• The styling of your question can be improved. Please take the time to read the markdown help page. You'll better communicate your problem when you use the right formatting.
– user9660
Nov 6, 2015 at 16:21
• You know Gamma is a Function, see Gamma, and your [ ] are possibly ill placed.
– user9660
Nov 6, 2015 at 16:26
• If the code is evaluated in order you are using % incorrectly. (first usage tries to get MX from the My,Mz solution for example) I'd suggest actually assiging the NDSolve results to symbols and using those instead of relying on output history (%) Nov 6, 2015 at 18:31

## Update

george2079 pointed out that using approximate numbers created a problem with the solution for high values of t. Replacing the approximate value for gamma with a rational number, the problem is stable at high values, including infinity.

γ = 267500000;
t1 = 1/10;
t2 = 1/100;
b0 = 2;
m0 = 1/10;
b1 = 1/10000000;


DSolve works with your equations.

The first differential equation yields a constant value for mx[t] of zero so I will skip it.

The second differential equation gives:

solYZ =
DSolve[{my'[t] == γ*mz[t]*b1 - my[t]/t2,
mz'[t] == -γ*my[t]*b1 + (m0 - mz[t])/t1, my[0] == 0,
mz[0] == m0}, {my[t], mz[t]}, t]

(* {{my[t] -> (1/11501679980)
107 E^(-55 t - (Sqrt[20951] t)/4 -
1/4 (-220 + Sqrt[20951]) t) (1676080 E^((Sqrt[20951] t)/2) -
838040 E^(1/4 (-220 + Sqrt[20951]) t) -
18649 Sqrt[20951] E^(1/4 (-220 + Sqrt[20951]) t) -
838040 E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t) +
18649 Sqrt[20951]
E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t)),
mz[t] -> (1/11501679980)
E^(-55 t - (Sqrt[20951] t)/4 -
1/4 (-220 + Sqrt[20951]) t) (670432000 E^((Sqrt[20951] t)/2) +
239867999 E^(1/4 (-220 + Sqrt[20951]) t) -
2518780 Sqrt[20951] E^(1/4 (-220 + Sqrt[20951]) t) +
239867999 E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t) +
2518780 Sqrt[20951]
E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t))}} *)


We can create functions out of the solutions in a number of ways. Below is one method:

my[t_] = solYZ[[1, 1, 2]]

(* (1/11501679980)107 E^(-55 t - (Sqrt[20951] t)/4 -
1/4 (-220 + Sqrt[20951]) t) (1676080 E^((Sqrt[20951] t)/2) -
838040 E^(1/4 (-220 + Sqrt[20951]) t) -
18649 Sqrt[20951] E^(1/4 (-220 + Sqrt[20951]) t) -
838040 E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t) +
18649 Sqrt[20951]
E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t)) *)


and for mz[t]

mz[t_] = solYZ[[1, 2, 2]]

(* (1/11501679980)E^(-55 t - (Sqrt[20951] t)/4 -
1/4 (-220 + Sqrt[20951]) t) (670432000 E^((Sqrt[20951] t)/2) +
239867999 E^(1/4 (-220 + Sqrt[20951]) t) -
2518780 Sqrt[20951] E^(1/4 (-220 + Sqrt[20951]) t) +
239867999 E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t) +
2518780 Sqrt[20951]
E^((Sqrt[20951] t)/2 + 1/4 (-220 + Sqrt[20951]) t)) *)


Now we can plot it:

Plot[{my[t], mz[t]}, {t, 0, 1},
PlotStyle -> {{Thick, Blue}, {Thick, Red}}]


It appears to quickly approach a constant. If you take the limits you will find:

Limit[my[t], t -> ∞]
(* 428/27449 *)

Limit[mz[t], t -> ∞]
(* 1600/27449 *)

• FWIW NDSolve does not diverge like that for large t. (perhaps a bifurcation?) Nov 6, 2015 at 18:45
• I'm pretty sure there is a numerical precision issue with evaluation the analytic solution. If you make all the parameters exact ( gamma=267500000) this goes away. Nov 6, 2015 at 19:06
• @george2079 You are correct. Answer updated. Nov 6, 2015 at 20:10

Not sure if this is what you are looking for but it's a start:

γ = 2.675*10^8;
T1 = 100*10^-3;
T2 = 10*10^-3;
B0 = 2;
M0 = 1*10^-1;
B1 = 1*10^-7;

nds1 = NDSolve[{Mx'[t] == -Mx[t]/T2, Mx[0] == 0}, Mx, {t, 0, 0.5}]

nds2 = NDSolve[{My'[t] == γ*Mz[t]*B1 - My[t]/T2,
Mz'[t] == -γ*My[t]*B1 + (M0 - Mz[t])/T1, My[0] == 0,
Mz[0] == M0}, {My, Mz}, {t, 0, 0.5}, MaxSteps -> Infinity]

Plot[Evaluate[Mx[t] /. nds1], {t, 0, 0.5},
FrameLabel -> {Style[t, 14], Style[Subscript[M, x], 12]},
PlotTheme -> "Scientific"]

Plot[Evaluate[{Mz[t], My[t]} /. nds2], {t, 0, 0.5},
PlotTheme -> "Scientific", PlotRange -> All,
FrameTicks -> {{Automatic, All}, {Automatic, None}}]