Clear["Global`*"];α=0.9;β=0.2;γ=1.2;
f[x_,y_,z_]:=z+(y-α)*xg[x_,y_,z_]:=1-β*y-x^2p[x_,y_,z_]:=-x-γ*z
Q[a_,b_,c_,h_,N_]:=(u[0]=a;v[0]=b;w[0]=c;Do[{u[n+1]=u[n]+h*f[u[n]+h/2*f[u[n],v[n],w[n]],v[n]+h/2*g[u[n],v[n],w[n]],w[n]+h/2*p[u[n],v[n],w[n]]],v[n+1]=v[n]+h*g[u[n]+h/2*f[u[n],v[n],w[n]],v[n]+h/2*g[u[n],v[n],w[n]],w[n]+h/2*p[u[n],v[n],w[n]]],w[n+1]=w[n]+h*p[u[n]+h/2*f[u[n],v[n],w[n]],v[n]+h/2*g[u[n],v[n],w[n]],w[n]+h/2*p[u[n],v[n],w[n]]]},{n,0,N}]);
Q[1,3,2,1,1000]
X=Interpolation[Table[{n,u[n]},{n,0,1000}]]
Y=Interpolation[Table[{n,v[n]},{n,0,1000}]]
Z=Interpolation[Table[{n,w[n]},{n,0,1000}]]
ParametricPlot3D[{X[t],Y[t],Z[t]},{t,0,1000}]