My code

function [ xs,ys,zs] = prog( ro,e0,alfa1,fi1,N,z )
emed=0; %energía del medio
%valores iniciales donde choca el rayo con la muestra
xi=0.5;
yi=0;
zi=0.5;
%vectores donde guardaremos todos los puntos
xs=xi;
ys=yi;
zs=zi;

%estos valores serán usados para que cada fotón comience con los mismos
%valores de alfa y fi, ya que los valores originales se pierden en el while
alfaprime=alfa1;
fiprime=fi1;

for i=1:N %número de fotones a simular
    a=2;

    alfa1=alfaprime;
    fi1=fiprime;

    while a>1
%valores iniciales donde choca el rayo con la muestra

m=9.10938356*(10^-31); %masa del electron
c=3*(10^8); %velocidad de la luz
gamma=e0/(m*c^2);

re=2.82*(10^-15);

dsigkn=(((re^2)/2)/(1+gamma*(1-cos(alfa1))^2))*(1+((cos(alfa1))^2)+gamma^2*(1-((cos(alfa1))^2))/(1+gamma*(1-cos(alfa1))));

sigkn=2*pi*(re^2)*(((1+gamma)/gamma^2)*((2*(1+gamma)/(1+2*gamma))-(log(1+2*gamma))/gamma)+((log(1+2*gamma))/(2*gamma))-((1+3*gamma)/(1+2*gamma)^2));

re=2.82*(10^-15);
fsc=(7.297352566)*(10^-3);
f=fsc*z;
fipe=4*pi*(re^2)*(f^4)*(z/gamma);
if f>=0 && f<=0.05
    sigpe=fipe*(1.008+(1.926/gamma)+(2.107/(gamma^2)))*1.079;
else if f>=0.05 && f<=0.1
        sigpe=fipe*(0.704+(1.647/gamma)+(1.592/(gamma^2)))*1.079;
    else if  f>=0.1 && f<=0.15
             sigpe=fipe*(0.604+(1.547/gamma)+(1.411/(gamma^2)))*1.092;
        else if f>=0.15 && f<=0.20
                sigpe=fipe*(0.522+(1.46/gamma)+(0.985/(gamma^2)))*1.092;
             else if f>=0.20 && f<=0.25
                sigpe=fipe*(0.455+(1.392/gamma)+(1.102/(gamma^2)))*1.102;
                 end
            end
        end
    end
end

sig=(z*sigkn)+sigpe;
c=z*sigkn/(sig);

disp('Ccccccc')
disp(c)

r=rand;
s=-log(r)/(sig*ro);

disp('S')
disp(s)
disp('Alfa1')
disp(alfa1)
disp('Fi1')
disp(fi1)

b=2;
criterio=0;
x1=s*cos(fi1)*sin(alfa1)+xi;
y1=s*sin(fi1)*sin(alfa1)+yi;
z1=s*cos(alfa1)+zi;
xi=x1;
yi=y1;
zi=z1;
if (x1>1)||(y1>1)||(z1>1)
disp('nuevo fotón porque el otro se salió')
break
end
if (x1<0)||(y1<0)||(z1<0)
disp('nuevo fotón porque el otro se salió')
    break
end
r2=rand;
%se decide que clase de efecto ocurre y los valores correspondientes

if (0<r2)&&(r2<c)
                r=rand;
                fi=2*pi*r;
                ry=20;
                p=0;
                   while ry>p
                       rx=pi*rand;
                       p=dsigkn*sin(rx)*2*pi;
                       ry=7*(10^-29)*rand;
                   end
                   alfa=rx;
else
   emed=e0+emed;
   b=0; %se sale del while
end

a=b;
if a==0
disp('nuevo fotón porque el anterior fue absorbido')
break
end
%lo que sigue es evaluar los ángulos de la nueva trayectoria si es que
%ocurrió compton

%con esto se calculan los ángulos de rebote, que usaremos para repetir el
%proceso
alfa2=acos(cos(alfa1)*cos(alfa)+sin(alfa1)*sin(alfa)*cos(fi));
disp('Alfa2')
disp(alfa2)
fi2=fi1+asin((sin(alfa)*sin(fi))/sin(alfa2));
disp('Fi2')
disp(fi2)
alfa1=alfa2;
fi1=fi2;
disp('Siguiente Alfa1')
disp(alfa1)
disp('Siguiente Fi1')
disp(fi1)

ef=e0/(1+gamma*(1-cos(alfa1))); %energ{ia del fotón reflejado
tmax=e0*(2*gamma/(1+2*gamma)); %máxima energía transferida al electrón
emed=tmax-ef+emed;
e0=ef;
     end
end

scatter(xs,ys)

end