@mmcdara When i run the your code, i got some errors in the code, is the code ok for you? If it is ok, then i explain the problem more.
As you're saying, this kind of differential integral equations are very complicated and i appreciate your efforts for solving this. I have some suggestions and explanations to improve this code,i hope we altogether with the help of you and other experts. could solve the problem.
1-# for n > 0 the denominator of op S is equivalent to -3.947841762*10^5*n^2 = (628.3185308*abs(n))^2 and the
# numerator is equivalent to 628.3185308*abs(n)*exp(-628.3185308*abs(n)*tau):
# Thus, assuming n > 0:) but in this work for getting the following equation
for calculating k1, we used a Taylor expansion which it has a sum on n=-infinity..infinity, and term of (nu[n]) is a frequency so it should be positive so we considered abs in k1 formula, so I think we should consider n<0 as n>=0 ( n=-infinity..infinity, ). what do you think? how we can extend this code to consider entire domain?
2- the other approximation is considered in
sum((0.1000000000e-1*exp(-0.1000000000e-1*tau)-628.3185308*abs(n)*exp(-628.3185308*abs(n)*tau))/(-394784.1762*n^2+0.1000000000e-3), n = -infinity .. infinity)
I think so, it is true only for very small amounts in comparison the others, what do you think?is it general?
3-in equation (7)
RHS := eval(-G*rho(tau)+1/2(G-F), rho(tau)=U);
could you please put a * between 1/2 and (G-F)? and i think it can affect the output.
4-the initial value of rho=0.5;
5-the maximum value of t you want to solve to:
almost after 200 (or more) it 'maybe' be constant, for example for T=100 the final constant is 0.4975 (steady state).
thank you for the time you spend on this,