Help with this code (fixed pivot technique)

Hi, Regards, i need help with this code, in order to get the molecular weight distribution of the polymer. This is using the fixed pivot technique.

here, i define the optimal stimated parameters:

clear all

clc

Param=[4.3276103255301335  0.10215543561231799  12.346130949956718  0.04565511507929096  0.0074744247798137625  0.0022142903501092513  1.6957452439712215  0.008166212885461337  0.028435528710812962  6.156650857937292  2.0566837008694536  0.9088238699999999  0.8669780915448081  8.072715437027492  0.010453310496367572  5.260491048  5.325728554979126  6.837140753063702  0.03947867100809781];

Optim=[1.4989542634814161  1.4335418277418217  1.5744495703184183  0.2559519340674219  0.7975221812123433  1.3319125502728477  1.1153538963112686];

Sens=[0.5  0.973342635  2  8.827704916  0.87519757  1.2];

% Parameters for Fed-Batch Models

um=0.54*Param(1)*0.6*0.6*0.5*Optim(1);

Ksr=0.15*Param(2); 

Sm=0.3*Param(3)*Optim(2);

Csx=1.7640*Param(4)*0.59*25*Optim(3);

Rcsx=0.022*Param(5)*Optim(4);

Csp=1.7460*Param(6);

k1=0.14*Param(7)*1.5*Optim(5);

k2=0.74*Param(8);

Cnx=0.336*Param(9);

KL=0.1*Param(10);

k3=93*Param(11);

k4=85*Param(12);

O2Leq=7.6*10^-3*Param(13);

nk=1.22*Param(14)*Optim(6);

Rcnx=0.16*Param(15)*Optim(7);

alfa1=0.143*Param(16);

alfa2=0.0000004*Param(17); 

alfa3=0.0001*Param(18); 

Kox=0.02*Param(19);

Now i define the state variables and feeding strategy for fed-batch fermentation solved by the Euler´s method

%Parameters for polymerization model

ki=0.62*10^4;

kp=0.46*10^5;

kt=0.14*10^1;

km1=0.11*10^-3;

km2=0.86*10^7;

kd=0.83*10^2;

%Initial Conditions validation data

X(1)=1.16;

S(1)=20.05;

N(1)=2;

P(1)=1.01;

O2L(1)=0.002432;

CO2(1)=0.01;

V(1)=4;

V1(1)=0;

V2(1)=0;

%Initial Conditions for polymerization

M(1)=0;

ESHM(1)=0;

dp1dt(1)=0;

dp2dt(1)=0;

dDdt(1)=0;

F3=0.1;

Sin=300*Sens(3);

Nin=7*Sens(4);

O2in=0.002432;

Yms=3.48*10^-3;

Jf=0.2818; 

Jm=Yms*Jf;

%Specific growth rate

u(1)=um*((N(1)/S(1))/((N(1)/S(1))+Ksr))*(1-((N(1)/S(1))/Sm)^nk)*((O2L(1)/(Kox*X(1)+O2L(1))));

tsim=50;

t(1)=0;

dt=0.001;

i=1;

Nt=2;

%Euler Method for ODE´s solution

while t(i)<tsim

      if t(i)<6
          F1=0;
          F2=0;
      elseif t(i)>6 && t(i)<10
          F1=0;
          F2=70*(1/1000)*Sens(5);
          Fobj31=1.3*F2*4;
      elseif t(i)>10 && t(i)<16
          F1=80*(1/1000)*Sens(6);
          F2=70*(1/1000)*Sens(5);
          Fobj32=3.2*F1*6+1.3*F2*6;
      elseif t(i)>16 && t(i)<20
          F1=80*(1/1000)*Sens(6);
          F2=0;
          Fobj33=3.2*F1*4;
          Fobj3=Fobj31+Fobj32+Fobj33;
      else
          F1=0;
          F2=0;

      end

      V(i+1)=V(i)+(F1+F2)*dt;
      V1(i+1)=F1;
      V2(i+1)=F2;

      %     Fed-Batch States
      X(i+1)=X(i)+(u(i)*X(i)*dt-((F1+F2)/V(i))*X(i)*dt);%Biomass
      S(i+1)=S(i)-((Csx*u(i)*X(i))+(Rcsx*X(i))+Csp*((k1*u(i)*X(i))+(k2*X(i))))*dt+(F1/V(i))*Sin*dt-((F1+F2)/V(i))*S(i)*dt;%Carbon Source
      P(i+1)=P(i)+((k1*u(i)*X(i))+(k2*X(i)))*dt-((F1+F2)/V(i))*P(i)*dt;%Polymer
      N(i+1)=N(i)-((Cnx*u(i)*X(i))+(Rcnx*X(i)))*dt+(F2/V(i))*Nin*dt-((F1+F2)/V(i))*N(i)*dt;%Nitrogen Source
      if N(i+1)<0.15
          N(i+1)=0.15;
      end
      O2L(i+1)=O2L(i)+((KL*(O2Leq-O2L(i)))-(k3*u(i)*X(i))-((k4*k1*u(i)*X(i))+(k4*k2*X(i))))*dt+(F3/V(i))*O2in*dt-((F1+F2)/V(i))*O2L(i)*dt;%Dissolved Oxygen
      if O2L (i+1)<0.002432
          O2L(i+1)=0.002432;
      end
      CO2(i+1)=CO2(i)+((alfa1*u(i)+alfa2)*X(i)*dt)+(alfa3*dt)-((F1+F2)/V(i))*CO2(i)*dt;%Dissolved CO2

      u(i+1)=um*((N(i+1)/S(i+1))/((N(i+1)/S(i+1))+Ksr))*(1-((N(i+1)/S(i+1))/Sm)^nk)*((O2L(i)/(Kox*X(i)+O2L(i))));
      if u(i+1)<0
          u(i+1)=0;
      end

Now, i am getting problem with the next part

%Monomer and Complex definition

      dp1dt(1)=1;%

      for n=1:Nt

      M(i+1)=M(i)+ (-((F1+F2)/V(i))*M(i)*dt+(-((Csx*u(i)*X(i))+(Rcsx*X(i))+Csp*((k1*u(i)*X(i))+(k2*X(i))))*Yms)-(km1*M(i)*dt)-(km2*M(i)*dp1dt(n))*dt); 

      ESHM(i+1)=ESHM(i)+(-((F1+F2)/V(i))*ESHM(i)*dt+(km1*M(i)*dt)-(ki*(ESHM(i)*dt))); 

for j=1:Nt

          for k=1:j

              dp1dt(j)=(-((F1+F2)/V(j))*dp1dt(j)+(ki(ESHM))-(km2*(dp1dt(j)*M))+kp*(dp2dt(k)*A(j,k)-kt*dp1dt(j)));

              dp2dt(j)=(-((F1+F2)/V(j))*dp2dt(j)+(km2*(dp1dt(j)*M))-(kp(dp2dt(j)))); 

              dDdt(j)=(-((F1+F2)/V(j))*dDdt(j)+(kt*dp1dt(j))-(kd*dDdt(j))+(kd*(dDdt(k)*B(j,k)))); 

          end

      end

      end

       t(i+1)=t(i)+dt;

      i=i+1;
  end

besides, A(j,k) and B(j,k) need to be solved by the fixed pivot technique, i´ll appreciate all your help.

[Walter Roberson]

"i am getting problem with the next part" is not very specific. Is the code meowing and demanding to be fed Purina Code Chow? Do vandals keep coming by and defacing the code? Is the code refusing to pay its taxes?