This is a quite general question, but I am working with a system of partial differential equations in two variables. There is one time direction t and one spatial direction z and the numerical method is formulated by stepping forward in time. The problem is that I obtain instabilities, either at the endpoints or sometimes even in the interior, depending on which method I am using. I have tried to approximate the spatial derivatives with finite differences of different orders, Runge Kutta, Euler backwards/forwards, pseudospectral methods based on Chebychev polynomials, I have tried to put the different functions on different grids, different methods to step forward in time, countless of ways to rewrite the equations or the order derivatives are evaluated. Even though there are some improvements, there is always some instability left and I am becoming desperate. Is there any expert out there who could just list all possible ideas one can try to make time evolution of PDEs stable? I dont think there is any point in writing out the equations here since they are quite long and complex and would just intimidate, BUT I can say that they are non-linear, involve both first and second order derivatives in both z and t (but they can be formulated in the standard form ∂_t f=… by suitable redefinitions) and the system includes five different functions that will constitute the solution.