How to declare j in runge kutta orde 4 function? Why my output doesnt exist?

clc;clear;

function RK4
%input konstanta
delta=50;
K1= 10^-4;
Ko=0.1;
n=3;
Oa=10;
Pa=100;
mu_1=10^-3;
mu_2=10^-3;
mu_3=10^-3;
mu_o=10^-4;
mu_p= 10^-5;
K2=5*10^-4;
K3=10^-4;
gamma=75;

%input initial condition
M1(1)=10;
M2(1)=0;
M3(1)=0;

%input for time
t(1)=0;
dt=0.01; %time interval
t=0:dt:100; %time span

%input empty array
T=zeros(length(t)+1,1); %empty array for t
M1=zeros(length(t)+1,1); %empty array for M1
M2=zeros(length(t)+1,1); %empty array for M2
M3=zeros(length(t)+1,1); %empty array for M3
O=zeros(length(t)+1,1);
P=zeros(length(t)+1,1);

fM1=@(M1,T) M1(j+1)+(dt*(delta*M1(j+1)*(1-(M1(j+1)/gamma))-2*K1*M1(j+1)*M1(j+1)-M1(j+1)*(K2.*M2(j+1))-((Oa-n)*K3*M1(j+1)*M3(j+1))-((Pa-Oa)*Ko*M1(j+1)*O(j+1))-(mu_1*M1(j+1))));

for j = 1:length((t)-1) 
    
    %T(j+1)=T(j)+dt;
    M1 =M1(j)+1./(1+exp(-T(j)));
    
    k1M1 = dt*fM1(M1(j),dt);
    k2M1 = dt*fM1(M1(j)+k1M1/2, dt+k1M1/2);
    k3M1 = dt*fM1(M1(j)+k2M1/2, dt+k2M1/2);
    K4M1 = dt*fM1(M1(j)+k3M1,dt+k3M1);
    M1(j+1) = M1(j)+1/6*(k1M1+2*k2M1+2*k3M1+k4M1);
    T(j+1)=T(j)+dt;
    
   
end

figure 
plot (T,M1,'r','Linewidth',3)
xlabel('time')
ylabel('M1')

end

[Torsten]

Please write down the differential equation you are trying to solve. Your rk4 code is completely weird, and I can't decipher your ODE from the code you posted.

[cindyawati cindyawati]

Thank you for your response @Torsten. I want to write this equation use RK4

[Torsten]

And where are the differential equations for M_2,...,M_n ?

Your array M must be of size (nt,n) where nt is the number of output times and n is the number of differential equations.

Similarly, your function fM1 must return a vector of length n for dM_1/dt,dM_2/dt,...,dM_n/dt.

[cindyawati cindyawati]

this is differential equation for M2 and Mn

[Torsten]

First write a function "fun" in which you define your system of differential equations:

function dy = fun(t,y)
  n = ...;
  M = y(1:n);
  O = y(n+1);
  P = y(n+2);
  dy(1) = ...;   % equations for M_1
  for i = 2:n
    dy(i) = ...;  % equation for M_i (2<=i<=n)
  end
  dy(n+1) = ...; % equation for O
  dy(n+2) = ...; % equation for P
end

Then we'll see how you should use it for RK4.

[cindyawati cindyawati]

thank you @Torsten I will try

[cindyawati cindyawati]

I get a new error @Torsten, in my function

[cindyawati cindyawati]

are you sure in O = y(n+1)? If I use n = 3 so in O= y(3+1). It is not equal to differensial of O

[Torsten]

are you sure in O = y(n+1)? If I use n = 3 so in O= y(3+1). It is not equal to differensial of O

I don't understand. You solve differential equations for M1,...Mn,O and P. You must define a vector y that comprises these unknown functions, and you save M(1) in y(1), M(2) in y(2), M(n) in y(n), O in y(n+1) and P in y(n+2).

If you use n = 3, you save M(1) in y(1), M(2) in y(2), M(3) in y(3), O in y(4)=y(n+1) and P in y(5)=y(n+2).

dy(1) = delta*M(1)*(1-M(1)/gamma)-2*K(1)M(1)^2-M(1)*sum(K(2:n-1).*M(2:n-1))-(Oa-n)*K(n)*M(1)*M(n)-(Pa-Oa)*Ko*M(1)*O-mu1*M(1)

Can you continue ?

And don't forget to supply the numerical values for all parameters used in "fun", i.e. delta, gamma, K(1)-K(n), Oa, Pa,Ko,...

[cindyawati cindyawati]

oke I will try again, sorry if this question bothers you

[Torsten]

If the question would bother me, I would not answer to it.

But I'm always surprised why people like you who have very little knowledge about MATLAB don't use MATLAB's onramp training course first before trying to solve rather complicated programming tasks:

https://matlabacademy.mathworks.com/details/matlab-onramp/gettingstarted