Array indices must be positive integers or logical values.

 f2=@(t,cSrum) (c1*exp(lambda(b-t))*(flm/fpf)*cSrlm(i))-(c1*exp(lambda(b-t))*flm*cSrlm(i)*18)-(c1*exp(lambda(b-t))*(1-flm)*c*25)+(c2*y(i)*exp(lambda*(b-t))*(1-fc)*cSrcc(i)*1);
    
    k9 =f2(t(i), cSrum(i));
    k10=f2(t(i)+0.5*h,cSrum(i)+0.5*k9*h);
    k11=f2(t(i)+0.5*h,cSrum(i)+0.5*k10*h);
    k12=f2(t(i)+h,cSrum(i)+k11*h);
    
    %update y;
    cSrum(i+1)=cSrum(i)+h/6*(k9+2*k10+2*k11+k12);

Error in

crustalgrowthradiogenicSr>@(t,cSrum)(c1*exp(lambda(b-t))*(flm/fpf)*cSrlm)-(c1*exp(lambda(b-t))*flm*cSrlm(i)*18)-(c1*exp(lambda(b-t))*(1-flm)*c*25)+(c2*y(i)*exp(lambda*(b-t))*(1-fc)*cSrcc(i)*1)

Error in crustalgrowthradiogenicSr (line 58)

k9 =f2(t(i), cSrum(i));

[James Tursa]

There are lots of places where this could be going wrong. Please show us the entire code so that we know what all of the variables involved are.

[Utpalendu Haldar]

%parameters
h=0.001;            %step size
tfinal=4.567;       %solve from t=0 to tfinal

%initial condition
t(1)=0;             %time
y(1)=0;             %mass of continental crust
cSrlm(1)=19.9;      %conc of Sr, BSE
cSrum(1)=19.9;      
cSrcc(1)=0;
b=4.567;
c2=0.15;
lambda=0.277;
m=2.2*10^22;
flm=0.1;
fc=0.6;
fpf=0.7*2.47*10^-3;


c1=c2*m*exp(c2/lambda*exp(lambda*b))/(exp(c2/lambda*exp(lambda*b))-exp(c2/lambda))

%RK4 loop
for i=1:ceil(tfinal/h)
    
    f=@(t,y) c1*exp((lambda*(b-t)))-c2*y*exp(lambda*(b-t));
    
    k1=f(t(i),y(i));
    k2=f(t(i)+0.5*h,y(i)+0.5*k1*h);
    k3=f(t(i)+0.5*h,y(i)+0.5*k2*h);
    k4=f(t(i)+h,y(i)+k3*h);
    %update y
    y(i+1)=y(i)+h/6*(k1+2*k2+2*k3+k4)
     
    fSrcc=@(t,cSrcc) c1*exp(lambda*(b-t))*(flm*cSrlm(i)*18+(1-flm)*cSrum(i)*25)-c2*y(i)*exp(lambda*(b-t))*cSrcc*1;
   
    k5=fSrcc(t(i),cSrcc(i))
    k6=fSrcc(t(i)+0.5*h,cSrcc(i)+0.5*k5*h);
    k7=fSrcc(t(i)+0.5*h,cSrcc(i)+0.5*k6*h);
    k8=fSrcc(t(i)+h,cSrcc(i)+k7*h);
    %update y;
    cSrcc(i+1)=cSrcc(i)+h/6*(k5+2*k6+2*k7+k8)
    
    fSrlm=@(t,cSrlm) c2*y(i)*exp(lambda*(b-t))*flm*cSrcc(i)*1-c1*exp(lambda*(b-t))*flm/fpf*cSrlm;
    
    k13=fSrlm(t(i),cSrlm(i));
    k14=fSrlm(t(i)+0.5*h,cSrlm(i)+0.5*k13*h);
    k15=fSrlm(t(i)+0.5*h,cSrlm(i)+0.5*k14*h);
    k16=fSrlm(t(i)+h,cSrlm(i)+k15*h);
    
    cSrlm(i+1)=cSrlm(i)+h/6*(k13+2*k14+2*k15+k16);
    
    f2=@(t,cSrum) (c1*exp(lambda(b-t))*(flm/fpf)*cSrlm(i))-(c1*exp(lambda(b-t))*flm*cSrlm(i)*18)-(c1*exp(lambda(b-t))*(1-flm)*c*25)+(c2*y(i)*exp(lambda*(b-t))*(1-fc)*cSrcc(i)*1);
    
    k9 =f2(t(i), cSrum(i));
    k10=f2(t(i)+0.5*h,cSrum(i)+0.5*k9*h);
    k11=f2(t(i)+0.5*h,cSrum(i)+0.5*k10*h);
    k12=f2(t(i)+h,cSrum(i)+k11*h);
    
   
    cSrum(i+1)=cSrum(i)+h/6*(k9+2*k10+2*k11+k12);
    
    %update t
    t(i+1)=t(i)+h;
    
    
end

 %plot result
 plot(t,cSrcc)
 hold on
 plot(t,cSrum)
 grid on
 
 xlabel('Time  (Gyr)')
 ylabel('Conc.')
 
 
 

this is the entire code

Answer 1

[Ameer Hamza]

Change the line

f2=@(t,cSrum) (c1*exp(lambda(b-t))*(flm/fpf)*cSrlm(i))-(c1*exp(lambda(b-t))*flm*cSrlm(i)*18)-(c1*exp(lambda(b-t))*(1-flm)*c*25)+(c2*y(i)*exp(lambda*(b-t))*(1-fc)*cSrcc(i)*1);

to

f2=@(t,cSrum) (c1*exp(lambda*(b-t))*(flm/fpf)*cSrlm(i))-(c1*exp(lambda*(b-t))*flm*cSrlm(i)*18)-(c1*exp(lambda*(b-t))*(1-flm)*c1*25)+(c2*y(i)*exp(lambda*(b-t))*(1-fc)*cSrcc(i)*1);)*cSrcc*1;
                           %^ missing                                 ^ missing                              ^ missing       ^ c is not defined        ^ missing