Regarding incorporating piecewise function in ode using octave code

Question

Differential equation where M_v is a piecewise function

Differential equation where M_v is a piecewise function

Piecewise function

I have coded the differential equation in Octave using the RUnge-Kutta 4 (RK4) method. However, the result is not as desired due to the piecewise function being coded incorrectly in my view. Could anyone help to code piecewise function correctly to solve differential equation using rk4 method?

  %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