Need help debugging so I can accurately portray graph of ODE system that is numerically solved

1 view (last 30 days)
So I'm able to produce the following graph but I need help debugging so I can figure out why it doesn't look the way it should like in the following paper. I haver a feeling that it is due to the value of h, the step size, it may be too big. Everything else seems ok, but I need to be sure
Here is my graph:
Here is graph in research paper(only part a):
Here is my code:
function cells
%constants
epsilon=0.0001; %detachment rate
M=30; %total number of ligands
N=4; %max binding number
Ns1=20; %number of nuclei for first set
Ns2=10; %number of nuclei for second set
%computing time variable
T=500000; %length of integration
tstart=10^-4;
tend=10^5;
t=tstart;
n=500; %total time steps
logstart=log10(tstart);
logend=log10(tend);
lpm=[logstart:(logend-logstart)/n:logend]; %log plot mesh
iplot=1; %counter
tarray=zeros(1,1);
tarray(iplot)=t;
h=max(10^(-5),0.1*min(1,t));
%Numerical scheme for sigma<1/2
c=[Ns1;0;0;0;0]; %initial conditions
carray=zeros(N+1,1); %removed fixed number of columns for matrix in order to amend column vectors in loop dependent on t
carray(:,iplot)=c; %adding initial condition to first slot of array
m=M-sum([1 2 3 4 5]'.*carray(:,iplot));
while t<=tend %loop for top graph implemented with Classical RK4
z=RK4(c,h,epsilon,m);
iplot=iplot+1;
carray=[carray z];
m=M-sum([1 2 3 4 5]'.*carray(:,iplot));
t=t+h;
tarray=[tarray t];
c=z;
h=max(10^(-5),0.1*min(1,t));
end
length(tarray)
length(carray(2,:))
logt=log10(tarray);
figure(1)
plot(logt,carray(2,:),'g')
hold on
plot(logt,carray(3,:),'y')
plot(logt,carray(4,:),'b')
plot(logt,carray(5,:),'k')
hold off
%Classical RK4
function W4=RK4(c,h,epsilon,m)
X1=c;
X2=c+h*(1/2)*rhs(X1,epsilon,m);
X3=c+h*(1/2)*rhs(X2,epsilon,m);
X4=c+h*rhs(X3,epsilon,m);
W4=c+h*((1/6)*rhs(X1,epsilon,m)+(1/3)*rhs(X2,epsilon,m)+(1/3)*rhs(X3,epsilon,m)+(1/6)*rhs(X4,epsilon,m));
%System of ODE's
function dcdt=rhs(c,epsilon,m)
dcdt=zeros(5,1);
dcdt(1)=-m*c(1)+epsilon*c(2);
dcdt(2)=-m*c(2)-epsilon*c(2)+m*c(1)+epsilon*c(3);
dcdt(3)=-m*c(3)-epsilon*c(3)+m*c(2)+epsilon*c(4);
dcdt(4)=-m*c(4)-epsilon*c(4)+m*c(3)+epsilon*c(5);
dcdt(5)=-epsilon*c(5)+m*c(4);

Answers (1)

Christopher Arreola
Christopher Arreola on 10 Dec 2021
I have fixed this, it was due to m

Categories

Find more on Loops and Conditional Statements in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!