Euler’s method for systems of ODEs in Mathematica
Math 330 
David McClendon

Command 1 (Defining code for Euler's method for systems of 2 ODEs) :

From lecture we described the formulas defining Euler’s method for first-order systems of 2 ODEs of the form y' =Φ(t,y).  To implement Euler's method for systems of 2 ODEs, first execute the following command once (this command will need to be executed once every time you restart Mathematica before you can apply it):

euler2D[{f_, g_},{t_,t0_,tn_},{x_,x0_},{y_,y0_},steps_]:=Block[{told=t0,xold = x0, yold=y0,thelist={{t0,{x0,y0}}},t,x,y,h},h=N[(tn-t0)/steps];
Do[tnew=told+h;
xnew = xold + h*(f/.{t->told,x->xold,y->yold});
ynew=yold+h*(g/.{t->told,x -> xold,y->yold});
thelist=Append[thelist,{tnew,{xnew,ynew}}];
told=tnew;
xold = xnew;
yold=ynew,{steps}];
Return[thelist];]

This command, in some sense, created a "black box" called euler2D which can henceforth be used to produce output obtained by Euler' s method.  After the above command is run once, you can then type in and execute commands like this:

Command 2 (Executing Euler's method for systems of 2 ODEs) :

euler2D[{y-x,x+y}, {t, 0, 3}, {x,2},{y, 0}, 3]

This command performs Euler’s method for the system of ODEs x’ = 2x + y, y’ = y - x^2 where the initial point is ((Subscript[t, 0]), (Subscript[y, 0])) = (1,(2.5,-2)); the procedure stops when t=3; there are 40 steps used.  This command will produce a (very long) list of points starting with (1, (2.5,-2)), (1.05, (2.65,-2.4125)), and ending with ((Subscript[t, 40]), (Subscript[y, 400])) = (3, (-1715.71, -50140.6).  This is the same as the list of points ((Subscript[t, n]), (Subscript[y, n])) that you would produce if you did Euler’s method by hand.  If you only want to see the last number in the list, tweak the command by typing

euler2D[{2x+y, y - x^2}, {t, 1, 3}, {x,2.5},{y, -2}, 40][[41]]

The number in the double brackets should be one more than the number of steps.

Command 3 (Sketching the points obtained from Euler's method for systems of 2 ODEs) :

If you only care about a picture of the points you’d get from Euler’s method, you can use the following command.  It simply plots all the points (Subscript[y, n]) that you would obtain by executing Euler’s method.  Note that the (Subscript[t, n]) are not indicated in this plot. Since there are so many points, it appears as though they make a curve (since the points are really close together):

ListPlot[Transpose[
euler2D[{y-x,x+y}, {t, 0, 3}, {x,2},{y, 0}, 3]][[2]],
 PlotRange -> {{-5,7},{-2,12}},
AspectRatio -> Automatic,
AxesLabel -> {"x","y"}]

Again, in this command the starting point is (1, (2.5,-2)), you stop at t = 3, and there are 40 steps used in the computation so there are 41 points plotted.   Here, the -6 and 6 are the smallest and largest values of x on the picture, and the -16 and 2 are the smallest and largest values of y on the picture.  Always keep the [[2]] at the end of the second line; this tells Mathematica to throw away the ts and only graph the xs and ys.

Command 4 (Defining code for Euler's method for systems of 3 ODEs) :

For systems of 3 ODEs, run this once every time you restart Mathematica to create the black box “euler3D”:

euler3D[{f_, g_,k_},{t_,t0_,tn_},{x_,x0_},{y_,y0_},{z_,z0_},steps_]:=Block[{told=t0,xold = x0, yold=y0,zold = z0, thelist={{t0,{x0,y0,z0}}},t,x,y,z,h},h=N[(tn-t0)/steps];
Do[tnew=told+h;
xnew = xold + h*(f/.{t->told,x->xold,y->yold, z-> zold});
ynew=yold+h*(g/.{t->told,x -> xold,y->yold, z-> zold});
znew=zold+h*(k/.{t->told,x -> xold,y->yold, z-> zold});
thelist=Append[thelist,{tnew,{xnew,ynew,znew}}];
told=tnew;
xold = xnew;
zold = znew;
yold=ynew,{steps}];
Return[thelist];]

Command 5 (Executing Euler’s method for systems of 3 ODEs) :

After having run Command 4 once, you can run this to produce the list of points for Euler’s method for a 3x3 system:

euler3D[{x - z, y + x, y - x}, {t, 0, 3}, {x,2},{y, -2}, {z,1}, 40]

The system of ODEs here is x’ = x-z; y’ = y+x; z’ = y-x;  the initial value is (0, (2,-2,1));  you use 40 steps and stop when t = 3.

To get only the last point, run this (the number in double brackets is one more than the number of steps):

euler3D[{x - z, y + x, y - x}, {t, 0, 3}, {x,2},{y, -2}, {z,1}, 40][[41]]

Command 6 (Sketching the points obtained from Euler’s method for systems of 3 ODEs) :

This command works similarly to Command 3; the initial value is the same as the one in Command 5.  The third part of the PlotRange command specifies the range of z values you see on the picture:

ListPointPlot3D[Transpose[
euler3D[{x - z, y + x, y - x}, {t, 0, 3}, {x,2},{y, -2}, {z,1}, 40]][[2]],
 PlotRange -> {{-6,6},{-16,2},{-10,4}},
AspectRatio -> Automatic,
AxesLabel -> {"x","y","z"}]