How to code impulsive differential equations in octave?

clc;
clear all;
close all;

% step sizes
t0 = 0;
tfinal = 5;
h = 0.01;
n = ceil((tfinal-t0)/h)+1;

%initial conditions
s(1) = 0.5; 
theta = 0.01;
N0 = 2;
Ta = 0.7;
b=N0*(Ta-theta)+theta;

x1(1) = -5;
x2(1) = 3;
t(1) =0;

% functions handles
Fx1 = @(t,x1,x2) -1/2*x1-x1.^3+x1.*x2.^2-((x1.^2+x2.^2).^1/4).*sign(x1);
Fx2 = @(t,x1,x2) -1/2*x2-x2.^3-x2.*x1.^2-((x1.^2+x2.^2).^1/4).*sign(x2);
for N = 1:n
  
if (mod(N,N0)==0)
  s(N+1) = s(N)+b;
else
  s(N+1) = theta;
end

  for i = 1:ceil(s(N+1))
 t(i+1) = t(i)+h;


##  if (mod(i,N0)==0)
##    h = b;
##    t(i+1) = t(i)+h;
##  else
##    h = theta;
##    t(i+1) = t(i)+h;
##  end
     %updates x1 and x2
 k1x1 = Fx1(t(i),     x1(i),       x2(i)        );
 k1x2 = Fx2(t(i),     x1(i),       x2(i)        );
 k2x1 = Fx1(t(i)+h/2, x1(i)+h/2*k1x1, x2(i)+h/2*k1x2);
 k2x2 = Fx2(t(i)+h/2, x1(i)+h/2*k1x1, x2(i)+h/2*k1x2);
 k3x1 = Fx1(t(i)+h/2, x1(i)+h/2*k2x1, x2(i)+h/2*k2x2);
 k3x2 = Fx2(t(i)+h/2, x1(i)+h/2*k2x1, x2(i)+h/2*k2x2);
 k4x1 = Fx1(t(i)+h,   x1(i)+h*k3x1, x2(i)+h*k3x2);
 k4x2 = Fx2(t(i)+h,   x1(i)+h*k3x1, x2(i)+h*k3x2);
 

if (i== ceil(s(N+1)))
 x1(i+1) = 0.5*x1(i)+h/6*(k1x1 +   2*k2x1 + 2*k3x1 +k4x1);
 x2(i+1) = 0.8*x2(i)+h/6*(k1x2 +   2*k2x2 + 2*k3x2 +k4x2);
else
 x1(i+1) = x1(i)+h/6*(k1x1 +   2*k2x1 + 2*k3x1 +k4x1);
 x2(i+1) = x2(i)+h/6*(k1x2 +   2*k2x2 + 2*k3x2 +k4x2);
end
end
N = N+1;
end

%plot
plot(t,x1);
hold on;
plot(t,x2);
##xlim([0, tfinal]);
##ylim([-6, 6]);
drawnow();