关于使用八度代码将分段功能合并到ODE中

Question

微分方程中m_v是一个分段函数

dindialial equication

” > dinactial方程src =“ https://i.sstatic.net/sfm55.png” alt =“ dindinial方程，其中m_v是一个零件函数”>

零件函数

我使用Runge-Kutta 4（RK4）方法对八度的微分方程进行了编码。但是，由于我认为分段函数编码错误，因此结果并没有达到期望。谁能帮助正确编码分段函数以使用RK4方法求解微分方程？

  %vIRt code of the equations given in the manuscript 

function x = pieceWise(s)
  x = zeros (size (s));
  ind = s > 0;
  x(ind) = s(ind);
endfunction
% Defines variables
global tau_s taur_a beta_r I_ext J_inter J_intra Jr_a  I_0 lp; 

tau_s  = 10; taur_a = 83;  beta_r = 0.0175; Jr_a = 172.97; J_inter= 81; J_intra = 32.4; I_0= 0.29; I_ext = 20;
% step sizes
t0 = 0;
tfinal = 200;
h = 0.05;
n = ceil((tfinal-t0)/h)+1;

%initial conditions
s_r(1) = 0;
s_p(1) = 0.3556;
a_p(1) = 6.151;
a_r(1) = 0;
t(1) =0;
 
% functions handles
S_r = @(t,s_r,s_p,a_r,a_p) -(s_r/tau_s)+ beta_r*pieceWise(I_ext-J_intra*s_r-J_inter*s_p-a_r-I_0);
A_r = @(t,s_r,s_p,a_r,a_p) (-a_r+Jr_a*beta_r*pieceWise(I_ext-J_intra*s_r-J_inter*s_p-a_r-I_0))/taur_a;
S_p = @(t,s_r,s_p,a_p,a_r) -(s_p/tau_s)+ beta_r*pieceWise(I_ext-J_inter*s_r-J_intra*s_p-a_p-I_0);
A_p = @(t,s_r,s_p,a_p,a_r) (-a_p+Jr_a*beta_r*pieceWise((I_ext-J_inter*s_r-J_intra*s_p-a_p-I_0)))/taur_a;


for i = 1:n
  
  %updates time 
  t(i+1) = t(i)+h;
  
  %updates S_r, A_r, S_p and A_p
 k1S_r = S_r(t(i),     s_r(i),         s_p(i),   a_r(i), a_p(i));
 k1A_r = A_r(t(i),     s_r(i),         s_p(i),   a_r(i), a_p(i));
 k1S_p = S_p(t(i),     s_r(i),         s_p(i),   a_r(i), a_p(i));
 k1A_p = A_p(t(i),     s_r(i),         s_p(i),   a_r(i), a_p(i));
 
 k2S_r = S_r(t(i)+h/2, s_r(i)+h/2*k1S_r, s_p(i)+h/2*k1S_p, a_r(i)+h/2*k1A_r, a_p(i)+h/2*k1A_p);
 k2A_r = A_r(t(i)+h/2, s_r(i)+h/2*k1S_r, s_p(i)+h/2*k1S_p, a_r(i)+h/2*k1A_r, a_p(i)+h/2*k1A_p);
 k2S_p = S_p(t(i)+h/2, s_r(i)+h/2*k1S_r, s_p(i)+h/2*k1S_p, a_r(i)+h/2*k1A_r, a_p(i)+h/2*k1A_p);
 k2A_p = A_p(t(i)+h/2, s_r(i)+h/2*k1S_r, s_p(i)+h/2*k1S_p, a_r(i)+h/2*k1A_r, a_p(i)+h/2*k1A_p);
 
 k3S_r = S_r(t(i)+h/2, s_r(i)+h/2*k2S_r, s_p(i)+h/2*k2S_p, a_r(i)+h/2*k2A_r, a_p(i)+h/2*k2A_p);
 k3A_r = A_r(t(i)+h/2, s_r(i)+h/2*k2S_r, s_p(i)+h/2*k2S_p, a_r(i)+h/2*k2A_r, a_p(i)+h/2*k2A_p);
 k3S_p = S_p(t(i)+h/2, s_r(i)+h/2*k2S_r, s_p(i)+h/2*k2S_p, a_r(i)+h/2*k2A_r, a_p(i)+h/2*k2A_p);
 k3A_p = A_p(t(i)+h/2, s_r(i)+h/2*k2S_r, s_p(i)+h/2*k2S_p, a_r(i)+h/2*k2A_r, a_p(i)+h/2*k2A_p);
 
 k4S_r = S_r(t(i)+h, s_r(i)+h*k3S_r, s_p(i)+h*k3S_p, a_r(i)+h*k3A_r, a_p(i)+h*k3A_p);
 k4A_r = A_r(t(i)+h, s_r(i)+h*k3S_r, s_p(i)+h*k3S_p, a_r(i)+h*k3A_r, a_p(i)+h*k3A_p);
 k4S_p = S_p(t(i)+h, s_r(i)+h*k3S_r, s_p(i)+h*k3S_p, a_r(i)+h*k3A_r, a_p(i)+h*k3A_p);
 k4A_p = A_p(t(i)+h, s_r(i)+h*k3S_r, s_p(i)+h*k3S_p, a_r(i)+h*k3A_r, a_p(i)+h*k3A_p);
 
 s_r(i+1) = s_r(i)+h/6*(k1S_r +   2*k2S_r + 2*k3S_r +k4S_r);
 s_p(i+1) = s_p(i)+h/6*(k1S_p +   2*k2S_p + 2*k3S_p +k4S_p);
 a_r(i+1) = a_r(i)+h/6*(k1A_r +   2*k2A_r + 2*k3A_r +k4A_r);
 a_p(i+1) = a_p(i)+h/6*(k1A_p +   2*k2A_p + 2*k3A_p +k4A_p);
end

M_r = beta_r*pieceWise((I_ext-J_intra*s_r-J_inter*s_p-a_r-I_0));
M_p = beta_r*pieceWise((I_ext-J_inter*s_r-J_intra*s_p-a_p-I_0));

%plots
##plot(t,M_p)
##hold on
##plot(t,M_r)
##plot(t,s_r)
##hold on
##plot(t,s_p)
##hold on
##plot(t,a_p)
##hold on
##plot(t,a_r)
##title('a_p and a_r')
##legend('s_r', 's_p', 'position', 'best')
##xlim([0,200])

drawnow;

%print -deps rk4.eps