Let's suppose that y is 1x5 and that the parameter is A, as in your diagram. Let's also assume that the equaitons of the systyem are too complicated to solve analytically. Therefore we will have to rely on numerical simulation. Make vector A with size Nx1, containing all the A values you want to have on your bifurcaiton diagram. Simulate the system with A=A(1). After a suitable simulation duraiton, evalute y(:) to see if it has reached a steady state (1). If y has reached a steady state, save the steady state values of y(1:5) as one row of a Nx5 array, for example yss1() (for steady state solution 1). We will use array yss1 for plotting later. Repeat with all N values of A. If you find no steady state solution for a particular value of A, save a row of NaNs in yss1(), for that particular vlaue of A. After you have done all N values of A, plot yss1(:,1) versus A. If there are multiple steady state solutions , you will have to find them by repeating the process above, but with different initial conditions. Your knowledge of the system should help you in this regard. On this second pass, in which you start from different inital points than what you used before, save the steady state solutions in yss2(). When done, add the plot of yss2(:,1) to the plot.
- To see if the system has reached a steady state, check y(:,1), y(:,2), ..., y(:,5) for multiple time points at the end of the simulation, to see if the vector y has stabilized.
I am attaching a recent publication (here, but may be paywalled), in which we published experimental and theoretical bifurcation diagrams for a complex system. Figures 4-9 are bifurcation diagrams in which the parameter on the horizontal axis is the frequency of external stimulation. This was a study of the repsonse of a non-linear optical system to external laser modulation. We did not do the bifurcaiton analysis as described above, because "steady state" meant a phase-locked solution, which was not a steady state in the traditional sense.
