Too many errors in your code. I tried to clean some, but it still has trouble to figure out what y is:
function kk1
clear all
close all
options = odeset('RelTol',1e-9,'AbsTol', 1e-16);
a0=0.1;
b0=0.1;
p=0.2;
theta=30:70;
alpha=sqrt(1-0.25*p^2);
a=(alpha-sqrt(1-(p)^2))./sqrt(alpha^2-0.25*sin(theta));
sol =ode45(@q,[0 1200],[0.01 0.01 0.01 0.01 0.01 0], options)
x1=sol(1).y(1,:);
y1=sol(1).y(2,:);
plot(x,y)
end
% % output
% x1=sol(1).y(1,:);
% y1=sol(1).y(2,:);
% plot(x,y)
%
function dy=q(x, alpha)
f=sqrt(1+2*x);
b=(1+2*x)^-1.5;
theta=30:70;
p=0.2;
m=16*alpha.^2-0.50*sin(theta);
m=4*(2*alpha.^2-8*(alpha*sqrt(1-p^2)));
o=16*sqrt(alpha.^2-0.50*sin(theta));
A=y*(1+2.*x)^-1.5.*(2+a); %What is y here?
B=-1i*4*y*(0.633)^2*(alpha.^2-0.25*sin(theta))+2*sqrt(alpha^.2-0.25.*sin(theta))*(alpha-sqrt(1-(p).^2));
C=y*0.633^2*(m+n+o).*(alpha-2*sqrt(1-(p).^2));
D=- y*((2+a)^-1)*f^-2*(1.5+1i*(9/8)*f*b);
E=-0.25*0.633*0.633*y*(4-5*p^2+4*(1-exp(0.18/4)));
F=(1/8)*p^2*b0*a0*((1.70*sin(theta)+3.5*sqrt(alpha^2-0.25*sin(theta))+2*(alpha-sqrt(1-(p)^2))));
dy=A+B+C+D+E+F;
end
% end
