Discrepancy of numerical results from analytical solution for an nonlinear viscoelasticity model

I have written a Fortran code for nonlinear viscoelasticity according to the Hughes' recurrence algorithm based on a total Lagrangian formulation. And a relaxation test is made. There is some deviation from the analytical solution(Theoretically the stress response should be an exponential decay plus a constant, isn't it?). The numerical ultimate stress is larger than the analytical one. As a larger initial strain is applied, the deviation is smaller. For a inital strain of 10%, there is almost no difference between the num. and analy. solution. What can account for this?