```Date: Dec 4, 2012 9:13 PM
Author: Kitech
Subject: for loop ode45 with dynamic variables

Hi I am writing a small program that can create many different particles that move up in the y-direction spread out in the x-direction. The particles interact with each other with repelling forces and so the position of each particle with reference to each other is required in the ode45 calculation. I have looked over the code many many times but I am not sure where the fault is as the solution comes up as "NaN." I have not shown constants in the code:time_interval = 1000; %Time in microseconds    %The number of time intervals the particle movement is divided by    number_of_intervals = 10000/time_interval;        tstart=0;    tstop= time_interval*1e-6;    num_of_particles = 3;        w0{1} = [0;0;0;0];    sols{1} = transpose(w0{1});    x_adj_stable(1) = 0;    y_adj_stable(1) = 0;        break_vect(num_of_particles) = 0;    times{1} = 0;        %Initiate particle characteristics    for particle_num = 2:num_of_particles    w0{particle_num}  = [0;0;w0{particle_num - 1}(3) + 4*1e-6;0];    x_adj_stable(particle_num) = w0{particle_num}(3);    y_adj_stable(particle_num) = 0;    sols{particle_num} = transpose(w0{particle_num});    times{particle_num} = 0;    end    %Event options, Identifying the time intervals within the solution    options = odeset('Events', @events);        %The loop which connects the intervals of time    for n = 1 : number_of_intervals                for particle_num = 1:num_of_particles    %===========================================    [t_sol, y_sol] = ode45(@pm1dwoc, [tstart, tstop], w0{particle_num}, options);    %===========================================    %The number of steps within each interval    n_steps = length(t_sol);           for interval = 2: n_steps            times{particle_num} = [times{particle_num}; t_sol(interval)];            sols{particle_num} = [sols{particle_num}; y_sol(interval, :)];    end    %Setting the initial conditions of the next interval    w0{particle_num} = y_sol(n_steps,:);    y_adj(particle_num) = y_sol(n_steps, 1);    x_adj(particle_num) = y_sol(n_steps, 3);             end        tstart = tstart + time_interval*1e-6;    tstop = tstart + time_interval*1e-6;        for n = 1:num_of_particles    y_adj_stable(n) = y_adj(n);    x_adj_stable(n) = x_adj(n);    end    end%=============================================== function dwdt = pm1dwoc (t,w)        y=w(1);        vy=w(2);        x=w(3);        vx=w(4);                for particle_number = 1:num_of_particles          Dist(particle_number) = sqrt((x - x_adj_stable(particle_number))^2 + (y - y_adj_stable(particle_number))^2);                if (Dist(particle_number) > 0)            x_vect(particle_number) = (x - x_adj_stable(particle_number))/Dist(particle_number);            y_vect(particle_number) = (y - y_adj_stable(particle_number))/Dist(particle_number);        end        end                             dy_component = -qw*Vd/(m*G) + qw*(-qw)/(m*4*pi*epsilon0*(2*R+2*y)^2));         dx_component = 0;         for particle_number = 1:num_of_particles             if (Dist(particle_number) > 0)                dx_component = dx_component + ((qw^2)/(0.0001*m*4*pi*epsilon0*Dist(particle_number)^2))*x_vect(particle_number);                dy_component = dy_component + ((qw^2)/(0.0001*m*4*pi*epsilon0*Dist(particle_number)^2))*y_vect(particle_number);             end         end                  dwdt = [vy; dy_component; vx; dx_component];             end    function [eventvalue,stopthecalc,eventdirection] = events (t,w)        t= t*1e6;        end_time = 10000;                eventvalue = t - end_time/number_of_intervals;        stopthecalc = 1;        eventdirection = 1;                for num = 2:number_of_intervals            eventvalue = [eventvalue; t - num*(end_time/number_of_intervals)];            stopthecalc = [stopthecalc; 1];            eventdirection = [eventdirection; 1];        end            end end   Any help would be really appreciated!Thank YouSam
```