Consider the following equation in my model
b=r*b(-1)+d;
After solving the model (either the nonlinear or linearised version), I take r and d generated by dynare and manually calculate the series of b (call it b2) using the equation above (again, either the nonlinear or linearised version).
I find the difference between b2 and b increases quickly in the length of simulation and becomes nonnegligible when the length is greater than 2000.
Is there a way to reduce such errors?

Are you solving the problem at first order? The equation you posted is nonlinear. If r is bigger than 1, you have an unstable difference equation where approximation errors can accumulate and explode.

Thanks for your reply.
I tried solving both the nonlinear and linearised versions.
Indeed, given r>1, approximation errors ought to explode.
I was hoping that maybe there is a way to reduce the errors such that the accumulated errors are still small enough after maybe 3000 periods

Yes, if solving the whole model to first order, I calculate b2 using the linearised version of the above law of motion.
You can easily replicate my problem in many standard model by adding

a gov budget constraint where there is a lump-sum tax

I just realised that there is basically no solution to the problem I posted.
Given r>1, there is only one b(-1) such that b=rb(-1)+d is stationary.
as T->infty, arbitrarily small deviation from that b(-1) will make b(t) explosive. so there is no room for any kind of error, even when the approximation error in dyanre is as small as rb(-1)+d-b<1e-15

There are, however, many ways to walk around this problem. Depending on what we want to do, we can for example

use knowledge about the terminal condition and solve backwards.

solve forward for a small T, once b(T) equal approximately its steady state, take b(t>T)=b(steady-state)

make a guess about b(t) for every t, and fsolve such that r*b(-1)+d-b is approximately zero.