Solving ODE with badly conditioned matrices

I'm trying to solve:

D*dUdt = -K*U + F(t, U)

Where D and K are sparse matrices produced by FEM software. The problem is that:

condest(D)     =  Inf
condest(K)     =  2.5628e+37
mean(mean(D))  = -1.5572e-10
min(min(K))    = -1.1956e+07
max(max(K))    =  2.2808e+07
mean(mean(K))  = -6.9619e-09

So usual solvers such as ode15s or ode23tb didn't work. I tried several methods to improve the conditioning number, such as adding a small value to the diagonal, using incomplete LU factorization to create a preconditioner matrix, and then feeding the system to an ODE solver but to not avail.

Now I'm solving "manually" the ODE with a backwards Euler method, which means for every time step n I have to solve:

A*Un+1 = bn+1
With
A = D/dt + K
bn+1 = D/dt*Un + Fn+1

To solve the linear system, I've been experimenting with several sparse matrix solvers such as bicgstab, gmres ... Until now I haven't found one that gives a satisfying result

Any suggestions ?

I've attached matrices D and K and initial conditions IC if someone wants to give it a shot. Fn should be a vector with:

Fn = zeros(662, 1);
Fn(662) = penalty*( Un(662) - (1-exp(-2*10^-3*tn)) );
With penalty something like 10^9

U(1) should look something like this: </matlabcentral/answers/uploaded_files/33070/U1.png>

U(14) something like this: </matlabcentral/answers/uploaded_files/33071/U14.png>

Thanks in advance !