StyleList = {AbsoluteThickness[1.0],#}& /@ {Blue,Red,Green,Purple};
SetOptions[Plot, Frame->True, Axes->{True,False}, PlotStyle->StyleList, AspectRatio->1/GoldenRatio ];
SetOptions[LogPlot, Frame->True, Axes->{True,False}, PlotStyle->StyleList, AspectRatio->1/GoldenRatio ];
SetOptions[ParametricPlot, Frame->True, Axes->{True,False}, PlotStyle->StyleList, AspectRatio->1/GoldenRatio ];
SetOptions[ListPlot, Frame->True, Axes->{True,False}, PlotStyle->{{AbsolutePointSize[5],Blue},{AbsolutePointSize[5],Red}}, AspectRatio->1/GoldenRatio ];SetOptions[ListLogPlot, Frame->True, Axes->{True,False}, PlotStyle->{{AbsolutePointSize[5],Blue},{AbsolutePointSize[5],Red}}, AspectRatio->1/GoldenRatio ];

Let us consider this simple system of differential equations

x'(t)=-1000 x
y'(t)=1000 x - y

with initial conditions x(0)=1 and y(0)=1. The exact solution is

sol=DSolve[{x'[t]== -1000 x[t],y'[t]== 1000 x[t]-y[t],x[0]==1,y[0]==0},{x,y},t][[1]]

Let us solve the system using the Euler implicit method. Here we do 100 steps.

euler=With[{h=0.1/1000},
f[{t_,x_,y_}]:= {t+h,x-1000 x h,y + h(1000 x - y )};
NestList[f,{0,1,0},100]];

Show[ListPlot[{euler[[All,{1,2}]],euler[[All,{1,3}]]}],Plot[Evaluate[{x[t],y[t]} /. sol],{t,0,0.01}]]

Increasing a factor 10 the step size

euler=With[{h=1/1000},
f[{t_,x_,y_}]:= {t+h,x-1000 x h,y + h(1000 x - y )};
NestList[f,{0,1,0},10]];

Show[ListPlot[{euler[[All,{1,2}]],euler[[All,{1,3}]]}],Plot[Evaluate[{x[t],y[t]} /. sol],{t,0,0.01}]]

Increasing a factor 2 more the step size

euler=With[{h=2.1/1000},
f[{t_,x_,y_}]:= {t+h,x-1000 x h,y + h(1000 x - y )};
NestList[f,{0,1,0},10]];

Show[ListPlot[{euler[[All,{1,2}]],euler[[All,{1,3}]]}],Plot[Evaluate[{x[t],y[t]} /. sol],{t,0,0.02}]]

In an implicit scheme we need to solve the following equation at each time step.

Solve[{xn==x-1000 xn h,yn== y + h (1000 xn-yn)},{xn,yn}]

implicit=With[{h=10./1000},
g[{t_,x_,y_}]:= {t+h,(x/(1+1000 h)),-((-1000 h x-y-1000 h y)/((1+h) (1+1000 h)))};
NestList[g,{0,1,0},10]];

Show[ListPlot[{implicit[[All,{1,2}]],implicit[[All,{1,3}]]}],Plot[Evaluate[{x[t],y[t]} /. sol],{t,0,0.1}]]

Show[ListLogPlot[{implicit[[All,{1,2}]],implicit[[All,{1,3}]]},PlotRange->{10^-6,2}],LogPlot[Evaluate[{x[t],y[t]} /. sol],{t,0,0.1},PlotRange->{10^-6,2}]]

The algorithm does not reproduce correctly the transition between x and y but the asymptotic values are correctly reproduced and we can evolve the whole system on much longer timescales.

implicit=With[{h=1./1000},
g[{t_,x_,y_}]:= {t+h,(x/(1+1000 h)),-((-1000 h x-y-1000 h y)/((1+h) (1+1000 h)))};
NestList[g,{0,1,0},10]];

Show[ListPlot[{implicit[[All,{1,2}]],implicit[[All,{1,3}]]}],Plot[Evaluate[{x[t],y[t]} /. sol],{t,0,0.01}]]

Show[ListLogPlot[{implicit[[All,{1,2}]],implicit[[All,{1,3}]]},PlotRange->{10^-6,2}],LogPlot[Evaluate[{x[t],y[t]} /. sol],{t,0,0.01},PlotRange->{10^-6,2}]]

implicit=With[{h=0.5/1000},
g[{t_,x_,y_}]:= {t+h,(x/(1+1000 h)),-((-1000 h x-y-1000 h y)/((1+h) (1+1000 h)))};
NestList[g,{0,1,0},20]];

Show[ListPlot[{implicit[[All,{1,2}]],implicit[[All,{1,3}]]}],Plot[Evaluate[{x[t],y[t]} /. sol],{t,0,0.01}]]

Show[ListLogPlot[{implicit[[All,{1,2}]],implicit[[All,{1,3}]]},PlotRange->{10^-6,2}],LogPlot[Evaluate[{x[t],y[t]} /. sol],{t,0,0.01},PlotRange->{10^-6,2}]]