# Euler's method for a 2nd order ODE

This is my first post on this site. Also, I'm new to Mathematica.

I'm trying to solve my first problem with Mathematica. It's about solving a 2nd order differential equation. I dont have the explicit F[t] instead I have a list of values for {t, F[t]}

 {{0, 1.00799}, {0.1, 1.09268}, {0.2, 1.18921}, {0.3, 1.25086}, {0.4, 1.32473},
{0.5, 1.36879}, {0.6, 1.39813}, {0.7, 1.41114}, {0.8, 1.39531}, {0.9, 1.3986},
{1., 1.39468}}


The equation I want to solve is

n = 5
h = 0.1
i = 0
y[0] = 1
yy[0] = 1
zi[0] = Table[i, {i, {1, 1}}]
T1 = {{0, 1.00799}, {0.1, 1.09268}, {0.2, 1.18921}, {0.3,
1.25086}, {0.4, 1.32473}, {0.5, 1.36879}, {0.6, 1.39813}, {0.7,
1.41114}, {0.8, 1.39531}, {0.9, 1.3986}, {1., 1.39468}}
Tij = Table[T1[[i, 2]], {i, 1, 11}]
ij = Table[T1[[i, 1]], {i, 1, 11}]
F[t_] := {y[t], (Tij[[t + 1]] - 9*yy[t] - 13*y[t])/5}
While[i < n, zi[i + 1] = zi[i] + F[i]*h; y[i + 1] = h + i;
yy[i + 1] = zi[[2]]; i = i + 1]

• Try again to paste your code. Do this: Put all your working code in one cell. Then do Cell->ConvertTo->InputForm (from the menu), then select the cell again and do RIGHT-CLICK (with the mouse), COPY-AS->Plain text. Then paste the result here again. Jun 15, 2014 at 7:49
• thx for that :) Jun 15, 2014 at 7:53
• I can't follow your code now. But noticed you wrote zi[[2]] there, and before you wrote zi[0, 0] You can't do this. zi[[2]] is trying to access position 2 in a list. You never allocated a list of this size. zi[2] is different from z[[2]]. You can write z[1000] without having to allocate anything (this is called indexed variable), but you can't write z[[1000]] without first allocating the space, as this is an actual list or vector in Matlab talk. Jun 15, 2014 at 8:02
• I think thats the problem... how can i enter a indexed vector like in this case zi[0]={1,1}, the indexed seems that only admits a lonely value and not vectors. (like y[0]=1) Jun 15, 2014 at 15:23

The easiest way to work with your sampled points is to derive an interpolating function from it.

pts =
{{0, 1.00799}, {0.1, 1.09268}, {0.2, 1.18921}, {0.3, 1.25086}, {0.4, 1.32473},
{0.5, 1.36879}, {0.6, 1.39813}, {0.7, 1.41114}, {0.8, 1.39531}, {0.9, 1.3986},
{1., 1.39468}};
f = Interpolation[pts];


Given f, to solve your ODE with Euler's method implemented by a while-loop, you could use

yPts = {{0., 1.}}; v0 = 1.;
Module[{t = yPts[[1, 1]], dt = .05, a, v = v0, y = yPts[[1, 2]]},
While[t <= 1.,
a = (9/5) v - (13/5) y + f[t]/5;
t += dt;
v += a dt;
y += v dt;
yPts = {yPts, {t, y}}];
yPts = Partition[Flatten[yPts], 2]];


Here is a plot of the results

yPlot = ListLinePlot[yPts, PlotStyle -> Red]


Mathematica, of course, can solve your ODE directly with NDSolve. It would be good to compare the above results with its more exact solution.

sol =
NDSolve[{y''[t] == (9/5) y'[t] - (13/5) y[t] + f[t]/5, y[0] == 1, y'[0] == 1},
y[t], {t, 0, 1}];