You first have to do some math:
Setting each of the functions to zero gives you two functions y(x) (called the nullclines), which you can plot in a phase diagram. Where the two lines intersect are the fixed-points (equilibria) of your system.
Now, you have to take the jacobian of your system and plug each of those fixed-points in, which will give you the linear stability analysis of the system.
The location of the fixed points and the stability of each point can now be computed as a you vary r (the bifurcation parameter).
For the programming:
-use newton's method (fsolve in MATLAB) to find where the equations are zero
-eig will help you find the eigenvalues of the system.
However
It depends on your system.
If you're supposed to be looking for limit cycles or chaos or something, you'll have to use one of the ode solvers and then the analysis becomes more tricky. I suppose you could develop a poincare-bendixson algorithm, but that would be involved and details would depend on your system.