Why am I not getting desired output?

I am solving differential equation that contain piecewise function. I coded this in Octave and I am getting desired results but the same in MATLAB is not giving the desired results. Code is compiling without any error but output is different. Code is given below:

function neural_circuit
% using rk4 method 
clc;
clear all;
close all;




h = 0.05;
t = 1:h:500;



y0 = [2.5, 3.6, 9.5, 8.8, 2, 3, 7, 7.5];
% t(1) =0;
y = zeros(length(t),8);
y(1,:) = y0;
for i = 1:length(t)

%   %updates time
   t(i+1) = t(i)+h;

    k1 = model(t(i), y(i, :));
    k2 = model(t(i)+h/2, y(i, :)+k1*h/2);
    k3 = model(t(i)+h/2, y(i, :)+k2*h/2);
    k4 = model(t(i)+h, y(i, :)+k3*h);
    y(i+1, :) = y(i, :)+(k1/6+k2/3+k3/3+k4/6)*h;

end

plot(t,y(:,1)) %% plot of s_l1

end


function dy= model(t,y)
%  dy = zeros(1,8);

%parameters value 
tau_s  = 10; taur_a = 83;  beta_l = 0.0175; Jr_a =148.20; I_ext = 20; I_0=0.29; J_l1l2 = 45.5; J_l1b1= 0.001; J_l2b2 = 0.001; J_b1b2= 44.5;
J_l1l1 = 35.4; J_l2l2 = 35.4; J_b1b1 = 32.4; J_b2b2 = 32.4;



function y = pieceWise(t)
    y = ...
       (t             ) .* (t >= 0) + ...
       (0 ) .* (t <0);
end

% function f = pieceWise(s)
%  if s>=0
%      f = s;
%  else
%      f=0;
%  end
% end

dy(1) = -(y(1)/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0);
dy(2) =  -(y(2)/tau_s)+ beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0);
dy(3) = -(y(3)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0);
dy(4) = -(y(4)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0);
dy(5) = -y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0)./taur_a;
dy(6) = -y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0)./taur_a;
dy(7) = -y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0)./taur_a;
dy(8) = -y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0)./taur_a;
end

[Jan]

clear all; on top of a function is a complete waste of time and has no meaningful purpose. It deletes all loaded functions from the RAM such that realoading them from the slow disk needs time without any benefit. This is called "cargo cult programming".

"J_l1l1" is really hard to read.

Answer 1

[Torsten]

h = 0.05;
t = 1:h:500;
y0 = [2.5, 3.6, 9.5, 8.8, 2, 3, 7, 7.5];
% t(1) =0;
y = zeros(length(t),8);
y(1,:) = y0;
for i = 1:length(t)
  t(i+1) = t(i)+h;
  k1 = model(t(i), y(i, :));
  k2 = model(t(i)+h/2, y(i, :)+k1*h/2);
  k3 = model(t(i)+h/2, y(i, :)+k2*h/2);
  k4 = model(t(i)+h, y(i, :)+k3*h);
  y(i+1, :) = y(i, :)+(k1/6+k2/3+k3/3+k4/6)*h;
end
plot(t,y(:,1)) %% plot of s_l1

function dy = model(t,y)
  tau_s  = 10; taur_a = 83;  beta_l = 0.0175; Jr_a =148.20; I_ext = 20; I_0=0.29; J_l1l2 = 45.5; J_l1b1= 0.001; J_l2b2 = 0.001; J_b1b2= 44.5;
  J_l1l1 = 35.4; J_l2l2 = 35.4; J_b1b1 = 32.4; J_b2b2 = 32.4;
  dy(1) = -(y(1)/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0);
  dy(2) = -(y(2)/tau_s)+ beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0);
  dy(3) = -(y(3)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0);
  dy(4) = -(y(4)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0);
  dy(5) = (-y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0))./taur_a;
  dy(6) = (-y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0))./taur_a;
  dy(7) = (-y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0))./taur_a;
  dy(8) = (-y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0))./taur_a;
end
function y = pieceWise(t)
  y = (t             ) .* (t >= 0) + ...
      (0 ) .* (t <0);
end

[Rakesh]

Thank you for your support!

Could you let know where did I do mistake?

[Torsten]

dy(5) = -y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0)./taur_a;
dy(6) = -y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0)./taur_a;
dy(7) = -y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0)./taur_a;
dy(8) = -y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0)./taur_a;

instead of

dy(5) = (-y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0))./taur_a;
dy(6) = (-y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0))./taur_a;
dy(7) = (-y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0))./taur_a;
dy(8) = (-y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0))./taur_a;

Answer 2

[Jan]

neural_circuit

function neural_circuit
h = 0.05;
t = 1:h:500;

y0 = [2.5, 3.6, 9.5, 8.8, 2, 3, 7, 7.5];
% t(1) =0;
y = zeros(length(t),8);
y(1,:) = y0;
for i = 1:length(t)
    t(i+1) = t(i)+h;
    k1 = model(t(i), y(i, :));
    k2 = model(t(i)+h/2, y(i, :)+k1*h/2);
    k3 = model(t(i)+h/2, y(i, :)+k2*h/2);
    k4 = model(t(i)+h, y(i, :)+k3*h);
    y(i+1, :) = y(i, :)+(k1/6+k2/3+k3/3+k4/6)*h;
end

plot(t,y(:,1)) %% plot of s_l1
end
function dy = model(t,y)
tau_s  = 10; taur_a = 83;  beta_l = 0.0175; Jr_a =148.20; I_ext = 20; I_0=0.29; J_l1l2 = 45.5; J_l1b1= 0.001; J_l2b2 = 0.001; J_b1b2= 44.5;
J_l1l1 = 35.4; J_l2l2 = 35.4; J_b1b1 = 32.4; J_b2b2 = 32.4;

function y = pieceWise(t)
    y = (t             ) .* (t >= 0) + ...
        (0 ) .* (t <0);
end
dy(1) = -(y(1)/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0);
dy(2) = -(y(2)/tau_s)+ beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0);
dy(3) = -(y(3)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0);
dy(4) = -(y(4)/tau_s)+ beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0);
dy(5) = -y(5)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(2)-J_l1b1*y(3)-J_l1l1*y(1)-y(5)-I_0)./taur_a;
dy(6) = -y(6)+Jr_a*beta_l.*pieceWise(I_ext-J_l1l2*y(1)-J_l2b2*y(4)-J_l2l2*y(2)-y(6)-I_0)./taur_a;
dy(7) = -y(7)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(4)-J_l1b1*y(1)-J_b1b1*y(3)-y(7)-I_0)./taur_a;
dy(8) = -y(8)+Jr_a*beta_l.*pieceWise(I_ext-J_b1b2*y(3)-J_l2b2*y(2)-J_b2b2*y(4)-y(8)-I_0)./taur_a;
end

Output with Octave:

I do not see the problem.

[Rakesh]

Output should be oscillatory. This is what I am getting in Octave.

[Jan]

Please copy and paste the code from my answer to your Octave and run it.

I'm convinced, that you use a different code - e.g. your diagram ends and 200.

[Rakesh]

You are right. I have coded in octave the same problem using function handle. Please see the octave code below---> this code is giving correct output but I have to code the same by spiliting the whole program into global and local functions that's why i have choosed MATLAB. But in matlab i am not getting the desired results...could you help to split the following program into main and local functions in octave?

% synapses to each sub-population.
graphics_toolkit('gnuplot')
clc;
clear all;
close all;


%defining piecewise functions
function x = pieceWise(s)
  x = zeros (size (s));
  ind = s > 0;
  x(ind) = s(ind);
endfunction


% Defines variables
global tau_s taur_a beta_l I_ext I_0 J_l1l2 J_l1b1 J_l1b2 J_l2b1 J_l2b2 J_b1b2 Jr_a M_l1 M_l2 M_b1 M_b2 J_l1l1 J_l2l2 J_b1b1 J_b2b2; 

tau_s  = 10; taur_a = 83;  beta_l = 0.0175; Jr_a =148.20;  I_0=0.29; J_l1l2 = 43.5; J_l1b1= 0; J_l2b2 = 0; J_b1b2= 44.5; 
J_l1l1 = 35.4; J_l2l2 = 35.4; J_b1b1 = 32.4; J_b2b2 = 32.4;
I_ext = 20;
% step sizes
##I_ext = 0:0.001: 20;
t0 = 0;
tfinal = 100;
h = 0.05;
n = ceil((tfinal-t0)/h)+1;

%initial conditions
s_l1(1) = 2.5;
s_l2(1) = 3.6;
s_b1(1) = 9.5;
s_b2(1) = 8.8;
a_l1(1) = 2;
a_l2(1) = 3;
a_b1(1) = 7;
a_b2(1) = 7.5;
t(1) =0;

% functions handles
S_l1 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) -(s_l1/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*s_l2-J_l1b1*s_b1-J_l1l1*s_l1-a_l1-I_0);
S_l2 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) -(s_l2/tau_s)+ beta_l*pieceWise(I_ext-J_l1l2*s_l1-J_l2b2*s_b2-J_l2l2*s_l2-a_l2-I_0);
S_b1 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) -(s_b1/tau_s)+ beta_l*pieceWise(I_ext-J_b1b2*s_b2-J_l1b1*s_l1-J_b1b1*s_b1-a_b1-I_0);
S_b2 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) -(s_b2/tau_s)+ beta_l*pieceWise(I_ext-J_b1b2*s_b1-J_l2b2*s_l2-J_b2b2*s_b2-a_b2-I_0);
A_l1 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) (-a_l1+Jr_a*beta_l*pieceWise(I_ext-J_l1l2*s_l2-J_l1b1*s_b1-J_l1l1*s_l1-a_l1-I_0))./taur_a;
A_l2 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) (-a_l2+Jr_a*beta_l*pieceWise(I_ext-J_l1l2*s_l1-J_l2b2*s_b2-J_l2l2*s_l2-a_l2-I_0))./taur_a;
A_b1 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) (-a_b1+Jr_a*beta_l*pieceWise(I_ext-J_b1b2*s_b2-J_l1b1*s_l1-J_b1b1*s_b1-a_b1-I_0))./taur_a;
A_b2 = @(t,s_l1,s_l2,s_b1,s_b2,a_l1,a_l2,a_b1,a_b2) (-a_b2+Jr_a*beta_l*pieceWise(I_ext-J_b1b2*s_b1-J_l2b2*s_l2-J_b2b2*s_b2-a_b2-I_0))./taur_a;



for i = 1:n
  
  %updates time 
  t(i+1) = t(i)+h;
  
%updates S_l1, S_l2, S_b1, S_b2, A_l1, A_l2, A_b1 and A_b2
 k1S_l1 = S_l1(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
 k1S_l2 = S_l2(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
 k1S_b1 = S_b1(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
 k1S_b2 = S_b2(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));

 k1A_l1 = A_l1(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
 k1A_l2 = A_l2(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
 k1A_b1 = A_b1(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
 k1A_b2 = A_b2(t(i), s_l1(i), s_l2(i), s_b1(i), s_b2(i), a_l1(i), a_l2(i), a_b1(i), a_b2(i));
 % k2th 
 k2S_l1 = S_l1(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
 k2S_l2 = S_l2(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
 k2S_b1 = S_b1(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
 k2S_b2 = S_b2(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
 
 k2A_l1 = A_l1(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
 k2A_l2 = A_l2(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
 k2A_b1 = A_b1(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
 k2A_b2 = A_b2(t(i)+h/2, s_l1(i)+h/2*k1S_l1, s_l2(i)+h/2*k1S_l2, s_b1(i)+h/2*k1S_b1, s_b2(i)+h/2*k1S_b2, a_l1(i)+h/2*k1A_l1, a_l2(i)+h/2*k1A_l2, a_b1(i)+h/2*k1A_b1, a_b2(i)+h/2*k1A_b2);
 % k3
 k3S_l1 = S_l1(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
 k3S_l2 = S_l2(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
 k3S_b1 = S_b1(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
 k3S_b2 = S_b2(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
 
 k3A_l1 = A_l1(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
 k3A_l2 = A_l2(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
 k3A_b1 = A_b1(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
 k3A_b2 = A_b2(t(i)+h/2, s_l1(i)+h/2*k2S_l1, s_l2(i)+h/2*k2S_l2, s_b1(i)+h/2*k2S_b1, s_b2(i)+h/2*k2S_b2, a_l1(i)+h/2*k2A_l1, a_l2(i)+h/2*k2A_l2, a_b1(i)+h/2*k2A_b1, a_b2(i)+h/2*k2A_b2);
 %k4
 k4S_l1 = S_l1(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
 k4S_l2 = S_l2(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
 k4S_b1 = S_b1(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
 k4S_b2 = S_b2(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
 
 k4A_l1 = A_l1(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
 k4A_l2 = A_l2(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
 k4A_b1 = A_b1(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
 k4A_b2 = A_b2(t(i)+h, s_l1(i)+h*k3S_l1, s_l2(i)+h*k3S_l2, s_b1(i)+h*k3S_b1, s_b2(i)+h*k3S_b2, a_l1(i)+h*k3A_l1, a_l2(i)+h*k3A_l2, a_b1(i)+h*k3A_b1, a_b2(i)+h*k3A_b2);
 %
 s_l1(i+1) = s_l1(i)+h/6*(k1S_l1 +   2*k2S_l1 + 2*k3S_l1 +k4S_l1);
 s_l2(i+1) = s_l2(i)+h/6*(k1S_l2 +   2*k2S_l2 + 2*k3S_l2 +k4S_l2);
 s_b1(i+1) = s_b1(i)+h/6*(k1S_b1 +   2*k2S_b1 + 2*k3S_b1 +k4S_b1);
 s_b2(i+1) = s_b2(i)+h/6*(k1S_b2 +   2*k2S_b2 + 2*k3S_b2 +k4S_b2);
 
 a_l1(i+1) = a_l1(i)+h/6*(k1A_l1 +   2*k2A_l1 + 2*k3A_l1 +k4A_l1);
 a_l2(i+1) = a_l2(i)+h/6*(k1A_l2 +   2*k2A_l2 + 2*k3A_l2 +k4A_l2);
 a_b1(i+1) = a_b1(i)+h/6*(k1A_b1 +   2*k2A_b1 + 2*k3A_b1 +k4A_b1);
 a_b2(i+1) = a_b2(i)+h/6*(k1A_b2 +   2*k2A_b2 + 2*k3A_b2 +k4A_b2);
end
%plots

##%average spike rates
M_l1 = beta_l*pieceWise(I_ext-J_l1l2*s_l2-J_l1b1*s_b1-a_l1-I_0);
M_l2 = beta_l*pieceWise(I_ext-J_l1l2*s_l1-J_l2b2*s_b2-a_l2-I_0);
M_b1 = beta_l*pieceWise(I_ext-J_l1b1*s_l1-J_b1b2*s_b2-a_b1-I_0);
M_b2 = beta_l*pieceWise(I_ext-J_l2b2*s_l2-J_b1b2*s_b1-a_b2-I_0);

figure(1);
freqz(M_l1, 1)

figure(2);
freqz(M_l2,1)
##figure(1);
##plot(t,M_l1)
##hold on 
##plot(t,M_l2)
##
##figure(2);
##plot(t,M_b1)
##hold on
##plot(t,M_b2)



##
##figure(1);
##plot(t,s_l1)
##hold on
##plot(t,s_l2)
##hold on
##figure(2);
##plot(t,s_b1)
##hold on
##plot(t,s_b2)

##plot(t,a_l1)
##hold on
##plot(t,a_l2)
##hold on
##plot(t,a_b1)
##hold on
##plot(t,a_b2)
##legend('s_b1', 's_b2', 'position', 'best')
##title('For B_1 and B_2 populations')
%legend('s_l1', 's_l2', 's_b1', 's_b2', 'position', 'best')
##xlim([0,50])

drawnow();