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

1 view (last 30 days)
Christopher Arreola on 9 Dec 2021
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);

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