% Constant parameters
DeltaE = 0 ;
omega = 0.1 ;
kesi = 8e-5 ;
Omega = 1e-2 ;
q = 1e-5 ;
% Initial Values
z0 = 0.15 ;
phi0 = (9/5)*pi ;
% undriven
epsilon = 0 ;
%%%driven
%%%epsilon =1e-3 ;
fun=@(t,y)[omega*y(2) + y(2)/sqrt(1-y(2)^2)*q*cos(y(1)) + epsilon*cos(t);-sqrt(1-y(2)^2)*sin(y(1))+ kesi*(omega*y(2) + y(2)/sqrt(1-y(2)^2)*q*cos(y(1)) + epsilon*cos(t))];
tspan = [0 30]
y0 = [phi0;z0];
[t,y] = ode45(fun,tspan,y0);
plot(t,y(:,2))
Best wishes
Torsten.
