i want to solve this, where σ = 0.1, L0 = 0.5, k= 0.01, tc = 70E-9
and U(Φ0) is
so that the phase difference Φ0 vs t graph and potential graph U(Φ0) vs Φ0 look like this 
I m getting this
for Φ0 vs t my program is
=================================================================
ti = 0; % inital time
tf = 10E-5; % final time
tspan = [ti tf];
o =10E5;
k = 0.01; % critical coupling strength
L = 0.5;
s = 0.1;
tc = 70E-9; % photon life time in cavity
f = @(t,y) [(-(s^2)*k/tc)*sin(y - pi/2) + L*(s^2)/(2*tc)*sin(y + pi/2)/sqrt(1 + cos(y + pi/2))];
[T, Y] = ode45(f,tspan, 0);
Y = linspace(-3,3,length(Y));
U = zeros(length(Y),1) ;
for i = 1:length(Y)
U(i) = -(1E-6).*(Y(i)) - (2 + (0.1)^2 ).*((0.033)./(70E-9)).*cos(Y(i) - pi/2) + (0.5).*(((0.1)^2)./(70E-9)).*((1 + cos(Y(i) + pi/2))^(0.5));
end
plot(Y,U);
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
xlabel('\phi')
ylabel('U')
=======================================================================
for U( Φ0) vs Φ0
==========================================================================
ti = 0; % inital time
tf = 10E-5; % final time
tspan = [ti tf];
o =10E5;
k = 0.01; % critical coupling strength
L = 0.5;
s = 0.1;
tc = 70E-9; % photon life time in cavity
f = @(t,y) [(-(s^2)*k/tc)*sin(y - pi/2) + L*(s^2)/(2*tc)*sin(y + pi/2)/sqrt(1 + cos(y + pi/2))];
[T, Y] = ode45(f,tspan, 0);
Y = linspace(-3,3,length(Y));
U = zeros(length(Y),1) ;
for i = 1:length(Y)
U(i) = -(1E6).*(Y(i)) - (2 + (0.1)^2 ).*((0.033)./(70E-9)).*cos(Y(i) - pi/2) + (0.5).*(((0.1)^2)./(70E-9)).*((1 + cos(Y(i) + pi/2))^(0.5));
end
plot(Y,U);
ax = gca;
ax.XAxisLocation = 'origin';
ax.YAxisLocation = 'origin';
xlabel('\phi')
ylabel('U')
====================================================================================
please tell me what's wrong in this program and what parameter can make my result correct?
