p1:=-1.31*^-2;
p2:=-1.35*^-2;
p3:=2.9*^-6;
h:=136;
n:=.13;
a:=3.11;
beta:=1;
gamma:=5.36*^-6;
G0:=200;
X0:=0;
I10:=0;
I20:=0;
Gb:=100;
Ib:=10;
tmax:=1000;
func:={(p1-X)*G-p1*Gb,p2*X+p3*(I1-Ib),a*Max[0,I2]-n*I1,beta*gamma*(G-h)-n*I2};
yinit:={G0,X0,I10,I20};
y:={G,X,I1,I2};
step:=.1
NDvec:={G'[t]==(p1-X[t])*G[t]-p1*Gb,X'[t]==p2*X[t]+p3*(I1[t]-Ib),I1'[t]==a*Max[0,I2[t]]-n*I1[t],I2'[t]==beta*gamma*(G[t]-h)-n*I2[t],G[0]==G0,X[0]==X0,I1[0]==I10,I2[0]==I20};
RungeKutta[func_List,yinit_List,y_List,step_]:=Module[{k1,k2,k3,k4},k1=step N[func/.MapThread[Rule,{y,yinit}]];
k2=step N[func/.MapThread[Rule,{y,k1/2+yinit}]];
k3=step N[func/.MapThread[Rule,{y,k2/2+yinit}]];
k4=step N[func/.MapThread[Rule,{y,k3+yinit}]];
yinit+Total[{k1,2 k2,2 k3,k4}]/6]
solutions:=NestList[RungeKutta[func,yinit,y,step]&,N[yinit],Round[tmax/step]]