How to vectorize this code?

clear
clc
tic
%concentration
%[M]=x1 [C]=x2 [Dr]=x3 [Cx]=x4 [R^o]=x5 [Rm^o]=x6 [Pr]=x7
x=[100 1 1 0 0 0 0];
v=1e-18;  % control volume
NA=6.022e23; % Avogadro's Number
t=0;
nmax=2000;
%number of molecules
n=(NA*v).*x;
N=n;
%kp=k1  ktc0=k2  ktd0=k3  ktr=k4   ka=k5 kd=k5
k=[1516 3.469e8 0 0.22 0.45 1.1e7];
kmc=k./((NA*v));
kmc2=(2.*k)./((NA*v));
iter=1;
iter_t=0;
for n=1:nmax
    for t=t
        R(1)=kmc2(5)*N(3)*N(2);
        R(2)=kmc2(6)*N(5)*N(4);
        R(3)=kmc2(1)*N(5)*N(1);
        R(4)=kmc2(4)*N(5)*N(1);
        R(5)=kmc2(3)*N(5)*N(6)+kmc(4)*N(5)*N(6);
        Rt=sum(R);
        P=[R(1)/Rt  (R(1)+R(2))/Rt (R(1)+R(2)+R(3))/Rt (R(1)+R(2)+R(3)+R(4))/Rt (R(1)+R(2)+R(3)+R(4)+R(5))/Rt];
        Pt=sum(P);
        r1=rand;
        r2=rand;
        tau=(1/Rt)*log(1/r2);
        meu=[1:1:5];
        if  r1<Pt
            %r3=rand;
            if (0<r1)&(r1<=P(1))
                
                N(2)=N(2)-1;
                N(3)=N(3)-1;
                N(4)=N(4)+1;
                N(5)=N(5)+1;
                
                
                % r4=rand
            elseif (P(1)<r1)&(r1<=P(2))
                
                N(5)=N(5)-1;
                N(2)=N(2)+1;
                N(4)=N(4)-1;
                N(3)=N(3)+1;
                
                % r5=rand
            elseif (P(2)<r1)&(r1<=P(3))
                
                N(1)=N(1)-1;
                N(5)=N(5)-1;
                N(6)=N(6)+1;
                
                % r6=rand
            elseif (P(3)<r1)&(r1<=P(4))
                
                N(1)=N(1)-1;
                N(7)=N(7)+1;
                
                %r7=rand
            elseif (P(4)<r1)&(r1<=P(5))
                
                N(5)=N(5)-2;
                N(7)=N(7)+2;
            end
            NN(iter,:)=N;
            N(1,:)=NN(iter,:);
            %
            
        else
            "mistake"
        end
        % iter=iter+1;
        %NN(iter,:)=N;
        %plot(t_save,NN(:,2))
    end
    %y=N(2)./(NA*v)
    NN(n,:)=N;
    N(1,:)=NN(n,:);
    Rt_save(n,:)=Rt;
    
    t_save(n,:) =t;
    tau_save(n,:)=tau;
    y=(N(1)./NN(:,1));
    y=vpa(y);
    z=100*(NN(:,1)/N(1));
    figure (1)
    plot(t_save,y,'-s')
    xlabel('polymerization Time(s)');
    ylabel('(M0/M)')
    figure (2)
    plot(t_save,z,'-o')
    xlabel('polymerization Time(s)');
    ylabel('Conversion')
    % hold on
    t=t+tau;
    
end
toc

[Geoff Hayes]

Taha - which part of your code do you want to vectorize? At the end of each iteration of the outer for loop you update a couple of figures...so if you vectorize part of this code (the outer loop) you might then lose the ability to visualize the changes on each iteration (unless you code for this). As for you inner for loop

for t=t

why bother having this? t will never change...

I would look closely at your code - some of what you do on each iteration of the outer loop doesn't even depend upon the n variable..unless you are resetting the N matrix? Without comments, it is unclear what should be happening.

Also, on each iteration of the loop you are calling plot. This will create a new graphics object each time so you may want to just create the object once and then update that object on subsequent iterations.

[Taha Abbas Bin Rashid]

Thanks for your time Geoff

variable 'n' is used to input total time of polymerization whereas 't' is changing after every iteration. After every iteration matrix 'N' also changes.

Regarding your suggestion of plot, I have done this and it resulted in huge time saving but I still need to make my code faster.

[Geoff Hayes]

Try preallocating or pre-sizing those arrays that you update on each iteration of the loop (this would be NN and N for example). See preallocation for details.

[Taha Abbas Bin Rashid]

Thanks a lot. Geoff