Thank you Roland for this detailed answer!
First of all, some things in my thorn have been changed in the last weeks. Therefore I attached the updated files to this mail. I reformulated the coefficients in my initial data. Now, the code runs much better but is still diverging. However the divergence seems to come from the center of the star instead. Right now I'm using the thornburgnc coordinates and made the spherical part as large as possible as this is the currently most stable configuration for my code.
I attached also a plot of the values of my evolved variable for a point at r=2 for two different time steps, where you can see the behavior of the evolution.
To be honest I'm not 100% sure what you meant with "stable solutions when implemented eg in spherical symmetry?"
What I'm basically trying to do is to have a wave equation that approximates the radial oscillation equations for small amplitudes. The coefficients of this wave equation are filled with static values using the TOVSolver. Then I only evolve the the differential equations, using the MoL thorn. There is no evolution of the GRHydro quantities or the ADM quantities themselves. In principal, time integration of these equations should be possible.
I wonder if I'm using the MoL (or other) thorn(s) wrong for this purpose?
Thank you a lot!
Best regards,
Severin