Ok, lasciamo perdere la cancellazione delle righe...
se io gli faccio risolve le eq. differenziali per ogni i e glidico di salvare i dati in T lui mi fa una matrice che non contiene le soluzioni...
usando questo ciclo:
for i=1:n              
         K_ion=normrnd(Ein,Ein/10);
         Etot_ion(i)=K_ion+((m_ion*uma1)/uma);  % ho tutto in MeV
         betasquare_ion(i)=1-(((m_ion*uma1)/(uma*Etot_ion(i)))^2);
         v_ion(i)=sqrt(betasquare_ion(i)*c^2);
         v_z=v_ion;
         
         %calcolo la componente x              
         v_x= unifrnd(-0.785,0.785, n,1); %angolo compreso da -45 e +45
                              % returns an array of random numbers generated from the continuous uniform distributions 
                              % with lower and upper endpoints specified by A and B, respectively. 
                              % If A and B are arrays, R(i,j) is generated from the distribution specified
                              % by the corresponding elements of A and B. If either A or B is a scalar, it is expanded
                              % to the size of the other input. R = unifrnd(A,B,m,n,...) or R = unifrnd(A,B,[m,n,...]) 
                              %returns an m-by-n-by-... array. If A and B are scalars, all elements of R are generated from 
                              %the same distribution. If either A or B is an array, they must be m-by-n-by-... 
          
          %calcolo la componente x                                  
          v_y= unifrnd(-0.785,0.785, n,1);
         
          v0=[v_x, v_y, v_z]; %initial velocity
          r0=zeros(n,3);      %initial speed
          
          y_sorgente=[r0,v0];  %vector of initial condition
          
          tspan=[0,100];
                    
          O=[0 0 0]';
          %%%Risolvo il moto
          f = @(t,ys) [ys(4:6); O];
          % the expression [y(4:6); (q/m)*cross(y(4:6),B)]
          % combines two vectors of length 3 to make a vector of
          % length 6(as long as y is a column vector).
          
          options=odeset('RelTol',1e-7,'Events',@(t,ys)Event_Stop_Sorgente(t,ys,Pinhole1));   % opzione per risoluzione equazione differenziale:
          %set precisione calcolo---% Sintassi per utilizzare Event_Stop con Pinhole1 come variabile in input
          [t,ys] = ode23t(f,tspan,y_sorgente(i,:),options);
          
          T(:,:,i)=ys(:,:); %Matrice tridimensionale delle traiettore di tutte le particelle
lui fa una T come quella in figura
mentre le soluzioni che calcola non hanno tutti elementi nulli a parte l'ultimo.
Allegati:
 9724_48f8b31ca23721529291553a81a070d4.jpg
9724_48f8b31ca23721529291553a81a070d4.jpg