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

a1=-1.0; b1=1.0;
a2=-1.0; b2=-1.;

x0=0.0;y0=0.0;tend=3.0;

sol=NDSolve[{
x'[t]==Cos[t],
y'[t]==Cos[t]+Exp[t]*(y[t]-Sin[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++,
{Print[solution=Solve[{
x1[n+1]==x1[n]+tau*Sin[(n+1)*tau],
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]/.solution[[1,1]];
x2[n+1]=x2[n+1]/.solution[[1,2]];
}]

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

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