function trajectory(x1,y1,z1,r)% spherical coordinate
for i=1:100
R=2*r;
theta=2*pi*(rand())^.5*((-1)^randi(2));;
Phi= 2*pi*(rand())^.5*((-1)^randi(2));;
x2= R*sin(theta)*cos(Phi);
y2= R*sin(theta)*sin(Phi);
z2= R*cos(theta);
plot3([x1,x2],[y1,y2],[z1,z2]);
hold on 
xlabel('x-axis')
ylabel('y-axis')
zlabel('z-axis')
grid on
axis ([ -10 10 -10 10])
end 
end