MATLAB: Please help me put this into a proper for loop, using increments

forloop

Hi, I am not very good with loops, so I did it be hand for 4 steps and then put them into code, I need somehow to put rFinal into a for loop so goes through 1000 times. I tried to do that, but without success, the loop are very incorrect, but I'm not sure how to fix them Use for reference, the variable status in order: Initial, Next, changed, New, Final Also NOTE: Pi in this case is the velocity not the variable "pi" (had to use variable pi..)
NOTE: the part I need help with is at the bottom (*** after the space), right BEFORE THE FOR LOOP, the rest is just variables if you want to run it…. if true %% Defining the variables
    % the masses
    M1=606; % galactic mass unit, = 2.32e7 of solar masses % central mass constants 
    M2=3690; % disk constants
    M3=4615; % halo constants
    % the radii and distances
    rInitial=8.499; % kpc  
    rMin=3.51; % kpc 

    rMax=9.15; % kpc 
    a2=5.3178; % disk constants % scalar lengths

    a3=12; % halo constants % scalar lengths
    b1=0.3873; % central mass constants % scalar lengths
    b2=0.25; % disk constants % scalar lengths
    z=0;
    % other remining variables and constants
    PiInitial=6.827; % kpc/gtm
    dt=1e-3; %time step
    %%the order of solving the variables - done by hand
    DeltarInitial= PiInitial*dt
    rNext= rInitial + DeltarInitial
    % solving for rChanged:
    % partials of Fi1 % the acceleration components - with r initial
    DiyFi1OverRInitial= (M1*rInitial) *dt / ( (b1^2+rInitial^2+z^2)^(3/2) );
    % partials of Fi2 % the acceleration components - with r initial
    DiyFi2OverRInitial= (M2*rInitial) *dt / ( (rInitial^2+( (a2+sqrt(b2^2+z^2))^2) )^(3/2) );
    % partials of Fi3 % the acceleration components - with r initial
    DiyFi3OverRInitialTerm1= ( 1.02*M3*rInitial*( (sqrt(rInitial^2+z^2) *dt / a3)^2.04 ) ) / (a3*(rInitial^2+z^2)*( 1+( (sqrt(rInitial^2+z^2) / a3)^1.02 )^2 ));
    DiyFi3OverRInitialTerm2= ( M3*rInitial*( ( sqrt(rInitial^2+z^2) / a3 )^2.02 ) ) *dt / ( ( (rInitial^2+z^2)^(3/2) )*( 1+( ( sqrt(rInitial^2+z^2) / a3 )^1.02 ) ) );
    DiyFi3OverRInitialTerm3= ( -2.02*M3*rInitial*( ( sqrt(rInitial^2+z^2) / a3 )^1.01 ) ) *dt / ( a3*(rInitial^2+z^2)*( 1+( ( sqrt(rInitial^2+z^2) / a3 )^1.02 ) ) ); % a negative term


    DiyFi3OverRInitial= DiyFi3OverRInitialTerm1 + DiyFi3OverRInitialTerm2 + DiyFi3OverRInitialTerm3;
    % the acceraration component in the r direction


    arInitial= DiyFi1OverRInitial + DiyFi2OverRInitial + DiyFi3OverRInitial
    PiNext=PiInitial + arInitial*dt
    DeltarNext= PiNext*dt
    rChanged= rNext + DeltarNext
    % sloving for rNew:
    % partials of Fi1 % the acceleration components - with rNext
    DiyFi1OverRNext= (M1*rNext) *dt / ( (b1^2+rNext^2+z^2)^(3/2) );
    % partials of Fi2 % the acceleration components - with rNext

    DiyFi2OverRNext= (M2*rNext) *dt / ( (rNext^2+( (a2+sqrt(b2^2+z^2))^2) )^(3/2) );
    % partials of Fi3 % the acceleration components - with rNext

    DiyFi3OverRNextTerm1= ( 1.02*M3*rNext*( (sqrt(rNext^2+z^2) *dt / a3)^2.04 ) ) / (a3*(rNext^2+z^2)*( 1+( (sqrt(rNext^2+z^2) / a3)^1.02 )^2 ));
    DiyFi3OverRNextTerm2= ( M3*rNext*( ( sqrt(rNext^2+z^2) / a3 )^2.02 ) ) *dt / ( ( (rNext^2+z^2)^(3/2) )*( 1+( ( sqrt(rNext^2+z^2) / a3 )^1.02 ) ) );
    DiyFi3OverRNextTerm3= ( -2.02*M3*rNext*( ( sqrt(rNext^2+z^2) / a3 )^1.01 ) ) *dt / ( a3*(rNext^2+z^2)*( 1+( ( sqrt(rNext^2+z^2) / a3 )^1.02 ) ) ); % a negative term
    DiyFi3OverRNext= DiyFi3OverRNextTerm1 + DiyFi3OverRNextTerm2 + DiyFi3OverRNextTerm3;
    % the acceraration component in the r direction
    arNext= DiyFi1OverRNext + DiyFi2OverRNext + DiyFi3OverRNext
    PiChanged=PiNext + arNext*dt
    DeltarChanged= PiChanged*dt
    rNew= rChanged + DeltarChanged
    %


    %
    %
    % sloving for rFinal: - THE LOOP PART******************************************
    % partials of Fi1 % the acceleration components - with rChanged
    DiyFi1OverRChanged= (M1*rChanged) *dt / ( (b1^2+rChanged^2+z^2)^(3/2) );
    % partials of Fi2 % the acceleration components - with rNext
    DiyFi2OverRChanged= (M2*rChanged) *dt / ( (rChanged^2+( (a2+sqrt(b2^2+z^2))^2) )^(3/2) );
    % partials of Fi3 % the acceleration components - with rNext
    DiyFi3OverRChangedTerm1= ( 1.02*M3*rChanged*( (sqrt(rChanged^2+z^2) *dt / a3)^2.04 ) ) / (a3*(rChanged^2+z^2)*( 1+( (sqrt(rChanged^2+z^2) / a3)^1.02 )^2 ));
    DiyFi3OverRChangedTerm2= ( M3*rChanged*( ( sqrt(rChanged^2+z^2) / a3 )^2.02 ) ) *dt / ( ( (rChanged^2+z^2)^(3/2) )*( 1+( ( sqrt(rChanged^2+z^2) / a3 )^1.02 ) ) );
    DiyFi3OverRChangedTerm3= ( -2.02*M3*rChanged*( ( sqrt(rChanged^2+z^2) / a3 )^1.01 ) ) *dt / ( a3*(rChanged^2+z^2)*( 1+( ( sqrt(rChanged^2+z^2) / a3 )^1.02 ) ) ); % a negative term
    DiyFi3OverRChanged= DiyFi3OverRChangedTerm1 + DiyFi3OverRChangedTerm2 + DiyFi3OverRChangedTerm3;
    % the acceraration component in the r direction
    arChanged= DiyFi1OverRChanged + DiyFi2OverRChanged + DiyFi3OverRChanged
    PiNew=PiChanged + arChanged*dt
    DeltarNew= PiNew*dt
    rFinal= rNew + DeltarNew
    for i=1:5
        arChanged= DiyFi1OverRChanged + DiyFi2OverRChanged + DiyFi3OverRChanged
        arChanged= arChanged + 1
        for i=1:5
            PiNew(i)=PiChanged + arChanged*dt
            for k=1:5
                DeltarNew(k)= PiNew*dt
                for i=1:5
                    rFinal(i)= rNew + DeltarNew                    
                    end  
                end
            end
        end
    end
end

Best Answer

The part after *************** can be written in loop as follow:
DiyFi1OverRChanged= (M1 * rChanged) * dt / ( (b1^2+rChanged^2+z^2)^(3/2) );
DiyFi2OverRChanged= (M2*rChanged) *dt / ( (rChanged^2+( (a2+sqrt(b2^2+z^2))^2) )^(3/2) );
DiyFi3OverRChangedTerm1= ( 1.02*M3*rChanged*( (sqrt(rChanged^2+z^2) *dt / a3)^2.04 ) ) / (a3*(rChanged^2+z^2)*( 1+( (sqrt(rChanged^2+z^2) / a3)^1.02 )^2 ));
DiyFi3OverRChangedTerm2= ( M3*rChanged*( ( sqrt(rChanged^2+z^2) / a3 )^2.02 ) ) *dt / ( ( (rChanged^2+z^2)^(3/2) )*( 1+( ( sqrt(rChanged^2+z^2) / a3 )^1.02 ) ) );
DiyFi3OverRChangedTerm3= ( -2.02*M3*rChanged*( ( sqrt(rChanged^2+z^2) / a3 )^1.01 ) ) *dt / ( a3*(rChanged^2+z^2)*( 1+( ( sqrt(rChanged^2+z^2) / a3 )^1.02 ) ) ); % a negative term
DiyFi3OverRChanged= DiyFi3OverRChangedTerm1 + DiyFi3OverRChangedTerm2 + DiyFi3OverRChangedTerm3;
arChanged= DiyFi1OverRChanged + DiyFi2OverRChanged + DiyFi3OverRChanged;
PiNew(1) = PiChanged + arChanged * dt;
DeltarNew(1) = PiNew(1) * dt;
rFinal(1) = rNew + DeltarNew(1);
for i=1:5
    arChanged= arChanged + 1;
    PiNew(i + 1)= PiChanged + arChanged * dt;
    DeltarNew(i + 1)= PiNew(i + 1) * dt;
    rFinal(i + 1)= rNew + DeltarNew(i + 1);
end
It it what you need?
Related Question