An earlier post discusses Mathematica code for Euler’s method. Here is some updated code.

A couple preliminary notes:

• I like to clear all variables at the beginning.
• For systems it is fun to superimpose the plot on top of the stream plot of the corresponding vector field

Example 1
For the IVP , with we can use the code

Clear["Global`*"]
f[t_, y_] := t - y
deltat = 0.1;
eulerValues = 
RecurrenceTable[{t[k + 1] == t[k] + deltat,
   y[k + 1] == y[k] + deltat*f[t[k],y[k]], 
   t[0] == 0., 
   y[0] == 2.}, 
   {t, y}, {k, 0, 10}];
Grid[eulerValues]
eulerPlot=ListLinePlot[eulerValues, PlotMarkers -> Automatic]
realPlot = Plot[3 Exp[-t] + t - 1, 
	{t, 0, 1}, PlotStyle -> Red];
Show[eulerPlot, realPlot]

Example 2
Consider , ,
We use

Clear["Global`*"]
deltat = .05;
f[y_] := Exp[y];
values = RecurrenceTable[{
    t[k + 1] == t[k] + deltat,
    y[k + 1] == y[k] + deltat*f[y[k]],
    t[0] == 0, y[0] == 1
    }, {t, y}, {k, 0, 10}];
Grid[values]
ListLinePlot[values, PlotMarkers -> Automatic]

Example 3
Here we consider the system


with and initial conditions . The following code processes the first 100 time steps.

Clear["Global`*"]

f[x_, y_] := x (1 - x) - x*y
g[x_, y_] := -y + 2 x*y
streamPlot = StreamPlot[
	{f[x, y], g[x, y]}, 
	{x, -1, 2}, 
	{y, -1, 2}]
	
deltat = 0.1;
steps = 100;

eulerValues = RecurrenceTable[{
    t[0] == 0, x[0] == 1.5, y[0] == 1,
    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}];
eulerPlotValues = Table[{
    eulerValues[[k, 2]],
    eulerValues[[k, 3]]},
   {k, 1, steps + 1}];
eulerPlot = ListLinePlot[
  eulerPlotValues,
  PlotMarkers -> Automatic,
  PlotStyle -> Red]
  
Show[streamPlot, eulerPlot]