hello everyone i have a problem to plot 3D nonlinear ode system ,this is my code and i don't know what is the error

% Oscillator function SOL=fnpp_new_1(t,X) x=X(1);y=X(2);z=X(3); alpha=10; beta=1; c=0.5; e=0.5; w=2/3; f=34.5; p= -1*((alpha^2)/4)+beta+ f*sin(z); Q= -1*(alpha^2/4)- beta/6 -( f*sin(z)/6); SOL(1)=y; SOL(2)=(p*x)-(c*y)-(e*y^3)+Q*(x^3); SOL(3)=w; SOL=SOL'; end clc;close all; clear; global f; f_arr=linspace(2,4,250); L_arr=zeros(0,length(f_arr)); step_s=800; %change this to see if it turns out well. for i=1:length(f_arr) f=f_arr(i); [t,L]=ode45(@fnpp_new_1,[0 5000],[0,0.2,2]); %fprintf('%d\n',i); size_temp=size(L(end-600:step_s:end,1)); f_arr2(1:size_temp)=f_arr(i); L_arr(1:size_temp)=L(end-600:step_s:end,1); plot (f_arr2,L_arr,'.'); hold on; %L_arr(1,400*(i-1)+1:400)=L(end-399:end,1); %this will plot bifurcation at x %clear t L; end %plot (a_arr2,L_arr,'.');

[nahla mahmoud]

by the way this is equation of my system dx/dt= y dy/dt= p*x-c*y-(e*(y^3))+Q*x.^3); p= -1*((alpha^2)/4)+beta+ f*sin(w*t); Q= -1*(alpha^2/4)- beta/6 -( f*sin(w*t)/6); alpha=10; beta=1; c=0.5; e=0.5; w=2/3; f=10;