this is a two-dimensional first order linear differential equation, and its general solution (without the constraints) is easy to compute
if we write x for the row vector (y1 y2), and A for the matrix (a c // b d), then the equation becomes
x' = x * A, x(0) = (p, r)
whose solution is
(see any introductory book on systems of differential equations)
x(t) = x(0) * exp(A * t)
in order to account for the constraints, we have to show that that solution is bounded as desired - that imposes strong constraints in the matrix A, expressible in terms of the eigenvalues of A - in other words, there is not always a solution that satisfies the constraints
because the constraints you have chosen mean that x must be bounded, i think that it is pretty straight forward to prove that the solution is periodic, and that the eigenvalues of A must be a complex pair with magnitude less than 1 (i have assumed that your coefficients and variables are all real numbers)
there is a marvelous book of Witold Hurewicz that explains all this (Lectures on ordinary differential equations)
The boundedness on y may suggest a periodic solution as you rightly said, this is on theoretical ground. Lets consider flow of liquid between two or tanks connected in different configurations, the simplest example is two tanks connected. modelling this with systems of first order ODE result in the given problem. the height of water is expected between two fixed height by controlling the flow. There are laboratory procedure for achieving this behavior but I need theoretical background in order to fix this for a large scale problems.
Ok, then i have a question. Are you starting with these constraints to find out whether there is a solution, or are you trying to find out whether or not there are bounds on a given problem?
As long as the equations are linear, you can do this kind of analysis in any finite number of dimensions, considering only the eigenvalues of the constant coefficient matrix A (there are many sources for this style of analysis; i'm using "stability theory of differential equations" by richard bellman, first published in 1953 and still useful)
A necessary and sufficient condition for all solutions to go asymptotically to 0 is that the real parts of all eigenvalues are negative. A necessary condition for bounded solutions is that the real part of each eigenvalue of A is not positive (Theorem 7 on p.25 and the exercise that follows it)
So it sounds like your original 2-d problem needs a real coefficient matrix A with two distinct eigenvalues whose real parts are non-positive, and at least one real part is 0. However, with a real coefficient matrix, the two eigenvalues are always complex conjugates, so to get one real part 0, the other one must also be 0, and therefore the two eigenvalues are both multiples of the pure imaginary number i, something like +w*i and -w*i for some real w > 0
Then the solutions are easy to compute
Of course, everything changes if the coefficient matrix can vary in time