# MATLAB: Runga Kutta numerical integration

oderunga kutta

Hi
I am trying to use the Runga Kutta method to solve 3 differential equations:
_dg/dt = [-g + f(i_c – i -h)]/b, where f is a function f(x)=0 if x<0, f(x)=x if 0<=x<a, f(x)=a if vm<=x
_dh/dt= (f*g -h)/d
_di/dt=(e*g -i)/c
where i_c, b, d,f,e and c are constants.
But this not yielding the expected result (a rapid rise in g when i_c changes from 1.15 to 3.1 followed by a rapid exponential decrease, and then a slow exponentatial decrease is expected) Does anyone know what might be wrong?
%clearclc;                                         clear;%constantsa = 4;b = 0.2;c = 24;d=294;e=6;f=9.4;k=1.15;l=3.1;%function handlesF_g=@(g,i_c,h,i,t) ((i_c -h -i)<0)*(-g/b)+ ((i_c -h -i)>=0)*((i_c -h -i)<a)*((i_c -h -i - g)/b) + ((i_c -h -i)>=a)*((a -g)/b);F_h=@(h,g,t) (f*g -h)/d;F_i=@(i,g,t) (e*g -i)/c;%steph=0.01;t = 0:h:1000; %prealocate memoryg=zeros(1,length(t));h=zeros(1,length(t));i=zeros(1,length(t));i_c=zeros(1,length(t));k_1_g=zeros(1,length(t));k_2_g=zeros(3,length(t));k_3_g=zeros(3,length(t));k_4_g=zeros(3,length(t));k_1_h=zeros(1,length(t));k_2_h=zeros(2,length(t));k_3_h=zeros(2,length(t));k_4_h=zeros(2,length(t));k_1_i=zeros(1,length(t));k_2_i=zeros(2,length(t));k_3_i=zeros(2,length(t));k_4_i=zeros(2,length(t));%initial valuesg(1)=k/(1+e+f);i(1)=e*g(1);h(1)=f*g(1);i_c(1:3000)=l;i_c(3001:5000)=k;i_c(5001:length(t))=l;for idx=2:(length(t))      k_1_g(idx) = F_g(t(idx-1),g(idx-1),h(idx-1),idx(idx-1));%v      k_1_h(idx) = F_h(t(idx-1),g(idx-1),h(idx-1));%as      k_1_i(idx) = F_i(t(idx-1),g(idx-1),i(idx-1));%af      k_2_g(idx) = F_g(t(idx-1)+0.5*h,g(idx-1)+0.5*h*k_1_g(idx),h(idx-1)+0.5*h*k_1_h(idx),i(idx-1)+0.5*h*k_1_i(idx));      k_2_h(idx) = F_h(t(idx-1)+0.5*h,h(idx-1)+0.5*h*k_1_h(idx),g(idx-1)+0.5*h*k_1_g(idx));      k_2_i(idx) = F_i(t(idx-1)+0.5*h,idx(idx-1)+0.5*h*k_1_i(idx),g(idx-1)+0.5*h*k_1_g(idx));      k_3_g(idx) = F_g((t(idx-1)+0.5*h),(g(idx-1)+0.5*h*k_2_g(idx)),(h(idx-1)+0.5*h*k_2_h(idx)),(i(idx-1)+0.5*h*k_2_i(idx)));      k_3_h(idx) = F_h((t(idx-1)+0.5*h),(h(idx-1)+0.5*h*k_2_h(idx)),(g(idx-1)+0.5*h*k_2_g(idx)));      k_3_i(idx) = F_i((t(idx-1)+0.5*h),(i(idx)+0.5*h*k_2_h(idx)),(g(idx-1)+0.5*h*k_2_g(idx)));      k_4_g(idx) = F_g((t(idx-1)+h),(g(idx-1)+k_3_g(idx)*h), (h(idx-1)+k_3_h(idx)*h), (i(idx-1)+k_3_i(idx)*h));      k_4_h(idx) = F_h((t(idx-1)+h),(h(idx-1)+k_3_h(idx)*h), (g(idx-1)+k_3_g(idx)*h));       k_4_i(idx) = F_i((t(idx-1)+h),(i(idx-1)+k_3_h(idx)*h),(g(idx-1)+k_3_g(idx)*h));      g(idx) = g(idx-1) + (1/6)*(k_1_g(idx)+2*k_2_g(idx)+2*k_3_g(idx)+k_4_g(idx))*h;      h(idx) = h(idx-1) + (1/6)*(k_1_h(idx)+2*k_2_h(idx)+2*k_3_h(idx)+k_4_h(idx))*h;      i(idx) = i(idx-1) + (1/6)*(k_1_i(idx)+2*k_2_i(idx)+2*k_3_i(idx)+k_4_i(idx))*h;% main equation  end  plot(t,g)