Here is a description of how to implement Euler’s method for a system.

We focus on a simple example:


where
and

First, let’s put in the functions which define the system

f[x_, y_] := x + y;
g[x_, y_] := -x;

Next, we make a StreamPlot of the vectorfield

sp = StreamPlot[{f[x, y], g[x, y]}, {x, -5, 5}, {y, -5, 5}];

Now we implement Euler’s method with and time steps by setting setting up three recursion relations:



The code which does this is the following:

deltat = .05;
steps = 100;
evalues = RecurrenceTable[{
    t[0] == 0,
    x[0] == 0,
    y[0] == .5,
    t[k + 1] == t[k] + deltat,
    x[k + 1] == x[k] + deltat*f[x[k], y[k]],
    y[k + 1] == y[k] + deltat*g[x[k], y[k]]
    },
   {t, x, y},
   {k, 0, steps}];
Grid[evalues];

Notice that I’ve put the semicolon after the Grid[evalues];. This suppresses the output. You can remove the semicolon to look at the list of numbers, which is helpful for finding errors.

Since we are only interested in plotting on the axis, we make a second table consisting only of the and values:

pvalues = Table[{
    evalues[[k, 2]],
    evalues[[k, 3]]
    }, {k, 1, steps + 1}];
Grid[pvalues];

You can use the Grid[pvalues] command to make sure you have the same list of numbers.

Now we can plot the numerical result on top of the StreamPlot:

eplot = ListLinePlot[pvalues, PlotMarkers -> Automatic, 
   PlotStyle -> Red];
Show[sp, eplot]