Bifurcation diagram for a system of 3 non-linear ODE

Hi everybody. I have a little problem to plot a bifurcation diagram for a system of 3 ODE. This is my system. The parameter is "da"

function dAdt = SistD( t,A )
global da mua mub nua nub nua_s etab xia xib
dAdt = zeros(3,1);
e=1e-2;
la=2.01+71*da;
la=la*e;
lb=-1.85+76.4i+(66.7+9.35i)*da;
lb=lb*e; 
dAdt(1)=la*A(1) -mua*A(1)*abs(A(1))^2 -nua*A(1)*abs(A(2))^2 -nua_s*A(1)*abs(A(3))^2 -xia*conj(A(1))*conj(A(3))*A(2);
dAdt(2)=lb*A(2) -mub*A(2)*abs(A(2))^2 -nub*A(2)*abs(A(3))^2 -etab*A(2)*abs(A(1))^2  -xib*A(3)*(A(1))^2;
dAdt(3)=lb*A(3) -mub*A(3)*abs(A(3))^2 -nub*A(3)*abs(A(2))^2 -etab*A(3)*abs(A(1))^2  -xib*A(2)*(conj(A(1)))^2;
end

In order to obtain the bifurcation diagram I solve the system with ode45. In this case I use the sum of the modulus of the three variables "phi" as "global variable" to understand the system behavior. I use the "stationary" value of the solution(the value that the solution reach for t>>0),the last 19 values of the solution, to plot the diagram.

clear variables
close all
clc
global da mua mub nua nub nua_s etab xia xib

e=1e-2;
%%Coeff A

mua=3.11;
mua=mua*e;
nua=6.88-1.11i;
nua=nua*e;
nua_s=6.88+1.11i;
nua_s=nua_s*e;
xia=4.57;
xia=xia*e;
%%Coeff B
mub=2.42+0.0321i;
mub=mub*e;
nub=3.13-0.816i;
nub=nub*e;
etab=0.955-3.47i;
etab=etab*e;
xib=1.62-1.36i;
xib=xib*e;

%%Calc
Rec=121.1;
eps=1e-1;
TF=500;
tspan=[0 700];
PHI=zeros(32,20);
options = odeset('RelTol',1E-12,'AbsTol',1E-15,'MaxStep',0.01);
%%Bif
for Re=115:147

   A0=[3 3 3]';
   da=(1/Rec-1/Re)*1/eps^2;

[t,L]=ode45(@SistD,tspan,A0,options);
TF=max(size(L));
phi=abs(L(:,1))+abs(L(:,2))+abs(L(:,3));
 PHI(Re-114,:)=phi(TF-19:TF)';
end
plot(115:1:147,PHI)

But I do not have the good solution. My questions are 1) Is my method right for this kind of non-linear system ? 2) The fact that i use conj function does affect the effectiveness of ode45? Thank you for your time and help

[Delfim Joao]

Hey, where you bale to solve this?