I didn’t see anything in your code that differed in any way from what the paper described. I didn’t read the paper as exhaustively as you must have, but I did read through it looking for a ‘Materials and Methods’ section, and didn’t find even the briefest mention of one. As you have seen with ode45 and ode15s, different solvers can give different results. The paper didn’t specify what solver they used. (It may have been the relatively popular and free Berkeley Madonna from U.C. Berkeley, but that’s simply a guess on my part. That is a fixed-step solver as best I can determine, but I haven’t examined it in detail because I have the MATLAB solvers.)
As for your getting complex results for x, the only possibility I can think of that could cause that is for (x(1) + x(2)) to become negative at some point. (I caution against using ‘i’ and ‘j’ as loop indices, because MATLAB uses them for its imaginary operators, though I don’t believe that’s causing those problems.)
I suggest you also look up the papers in the references by Gaffney, and Piccart-Gebhart that may help elucidate the differences you found. You can also write the authors.
That a paper is published in a well-respected, refereed journal does not guarantee that it is possible to reproduce the results, especially if there is no mention of methods.
Best Answer