Hi. Are you required to use this numerical integration method exactly ? Where did this limitation came from? Should you consider embedded ODE solvers of MATLAB, or construct algrebraic equations on yourself?

Below you can find a simplest Euler method code in C-like syntax for Hindmarsch-Rose model. It is possible to create any RK method from here.

x[0]=x1[0]+h*(x[1]+a[0]*x[0]**2-x[0]**3-x[2]+i);

x[1]=x1[1]+h*(a[1]-a[2]*x[0]**2-x[1]);

x[2]=x1[2]+h*(a[3]*(a[4]*(x[0]+a[5])-x[2]));

Where x[3] is the state variable array, x1[3] is the array of state variables values on the previous timestep, a=[3, 1, 5, 0.01, 4, 1.6] is the array of system parameters, and h is the integration stepsize.

In case you need to use embedded RK4(5) ODE solver in MATLAB, just rewrite the right-hand side without "x1[n]" and "h*" members and apply the default MATLAB solver.