a1=-500000.5; b1=499999.5;
a2=499999.5; b2=-500000.5;

a1=-1.0; b1=0.0;
a2=0.0; b2=-10.0;

x0=0.0;y0=2.0;tend=4.0;

sol=NDSolve[{
x'[t]==a1*x[t]+b1*y[t],
y'[t]==a2*x[t]+b2*y[t],
x[0]==x0,y[0]==y0},{x,y},{t,0,tend}]

Plot[Evaluate[{x[t],y[t]}/.sol],{t,0,tend}]

Clear[x1];Clear[x2];

x1[0]=x0; x2[0]=y0; nn=10;tau=tend/nn;

For[n=0, n<nn, n++,
{
sol=Solve
[
{
x1[n+1]==x1[n]+tau*(a1*x1[n+1]+b1*x2[n+1]),
x2[n+1]==x2[n]+tau*(a2*x1[n+1]+b2*x2[n+1])
},
{x1[n+1],x2[n+1]}
];
x1[n+1]=x1[n+1]/.sol[[1,1]];
x2[n+1]=x2[n+1]/.sol[[1,2]];
}]

ListLinePlot[Table[{n*tau,x1[n]},{n,0,nn}]]

ListLinePlot[Table[{n*tau,x2[n]},{n,0,nn}]]