# Euler's Method for system of differential equation with more than two functions

According to what I have read, the first formula is the classical Euler method and the second is the improved Euler method for second-order equations.

• Method A: accuracy of order h
S[a_, b_, h_, N_] := (u[0] = a; u[1] = a + h*b;
Do[u[n + 1] =
2 u[n] - u[n - 1] + h*h*f[n*h, u[n], (u[n] - u[n - 1])/h], {n, 1,
N}])

• Method B: accuracy of order h^2
Q[a_, b_, h_, N_] := (u[0] = a; v[0] = b;
Do[{u[n + 1] =
u[n] + h*
F[u[n] + (h/2)*F[u[n], v[n]],
v[n] + (h/2)*
G[u[n], v[
n]]],                                                     \

v[n + 1] =
v[n] + h*
G[u[n] + (h/2)*F[u[n], v[n]], v[n] + (h/2)*G[u[n], v[n]]]}, {n,
0, N}])


All the examples I have seen both in the book I have and in YouTube videos deal with two function systems and initial values. My question is this if we are given more than two functions how should we work? I can't find an example guide to figure out how. Any example would be appreciated. Thank you in advance

For example

It is given the following problem: $$X'=Z+(Y-\alpha)X$$ $$Y'=1-\beta Y-X^2$$ $$Z'=-X-\gamma Z$$ with initial conditions $$(X(0),Y(0),Z(0)=(1,2,3)$$.

Where

X: interest rate

Υ:investment demand

Z: price index

$$\alpha$$: savings, $$\beta$$: cost per investment, $$\gamma$$: the absolute value of the elasticity of demand

And we want the results with the Improved Eulers method. It is obvious that I have to use Method B: accuracy of order h^2. But I do not know how to define the new function on Mathematica

• To me, this sounds more of a Mathematics Question, appropriate for Math SE. Commented Jan 26, 2023 at 23:46
• @Domen Fyi I have updated my post so you will understand what I mean Commented Jan 26, 2023 at 23:56
• @cvgmt I created a new question about Improved Eulers Method as we discussed earlier today. Commented Jan 27, 2023 at 0:27
• 1. 1st implementation you've found is bounded on second order equations i.e. equation in form $y′′=f(x,y,y′)$ so it cannot be extended to first order system i.e. $y′=f(x,y)$, at least cannot be extended in a straightforward way. 2. Euler's method can be defined based on first order system in a much simpler manner, and n-th order system can always be transformed to 1st order system, see e.g. mathematica.stackexchange.com/a/158519/1871, notice 2nd implementation you've found is essentially based on this idea. BTW… Commented Jan 27, 2023 at 2:59
• …it's usually called midpoint method: en.wikipedia.org/wiki/Midpoint_method 3. 2nd implementation you've found isn't convenient for extending, to implement it in a more convenient way, observe the output of Sin[{Pi/3, Pi/6}] and think about how you can make use of this feature of Mathematica. 4. Please first make some effort to understand Euler's method itself, the wiki page isn't a bad source: en.wikipedia.org/wiki/Euler_method Commented Jan 27, 2023 at 3:07

This works but it is not as accurate as the NDSolve result

α=0.9;
β=0.2;
γ=1.2;

f[x_,y_,z_]:=z+(y-α)*x
g[x_,y_,z_]:=1-β*y-x^2
p[x_,y_,z_]:=-x-γ*z

Q[a_,b_,c_,h_,N_]:=(u[0]=a;v[0]=b;w[0]=c;
Do[{u[n+1]=u[n]+h*f[u[n]+h/2*f[u[n],v[n],w[n]],
v[n]+h/2*g[u[n],v[n],w[n]],
w[n]+h/2*p[u[n],v[n],w[n]]],
v[n+1]=v[n]+h*g[u[n]+h/2*f[u[n],v[n],w[n]],
v[n]+h/2*g[u[n],v[n],w[n]],
w[n]+h/2*p[u[n],v[n],w[n]]],
w[n+1]=w[n]+h*p[u[n]+h/2*f[u[n],v[n],w[n]],
v[n]+h/2*g[u[n],v[n],w[n]],
w[n]+h/2*p[u[n],v[n],w[n]]]},
{n,0,N}]);

Q[1,3,2,0.1,1000]

X=Interpolation[Table[{n,u[n]},{n,0,1000}]];

Y=Interpolation[Table[{n,v[n]},{n,0,1000}]];

Z=Interpolation[Table[{n,w[n]},{n,0,1000}]];

ParametricPlot3D[{X[t],Y[t],Z[t]},{t,0,1000}]


• Not too bad in my view. You can improve it by making the grid denser, try e.g. h = 0.01; tend = 100; nmax = tend/h; Q[1, 3, 2, h, nmax] // AbsoluteTiming Commented Jan 28, 2023 at 2:49