Solving differential equations on c++

The main function runs a loop that feeds vector values to a serie of sequentially interdependent differential equations found in another function gate_probabilities, which feeds back the calculations to main and the loop repeats.

The solutions for the equations are serially fed each into its corresponding vector, this takes place in both functions. intermediate equations are m_new; h_new; n_new, the main vector to eventually be plotted against time is V. (plotting not in the code)

The vectors are defined outside the function body (they are thought of as global but i guess the issue is they're not properly declared as such)

Building produces no errors, but i do get a debug error on execution: "Abort() has been called".

would appreciate any help :)

#include "C:\Users\jenbe\Documents\Visual Studio 2015\std_lib_facilities.h"
#include <math.h>
#include <iostream>
#include <cmath>
#include <tuple>
#include <vector>

using namespace std;

vector<double> mvector;
vector<double> hvector;
vector<double> nvector;
vector<double> V; // potential vector
vector<double> timevector;

//global constants

double gL = 0.1; // mS/cm^2
double gK = 9;
double gNa = 35;
double EL = -65; // mV
double EK = -90;
double ENa = 55;
double phi = 5; // Coefficient increasing reaction speed of the alpha and beta constants which it multiplies.
double Iapp = 5; // microAmpere values (as in Wang-Buzsaki)
double runtime = 3000;// 3 Seconds
double dt = 0.001; // Timestep = 1ms
double V_init = -60;


std::tuple<double> gate_probabilities(double v, double m, double h, double n) { 
    double am, bm, ah, bh, an, bn, dh, dn;

    //first stage -> gate probabilities

    am = -0.1*(v + 35) / (exp(-0.1*(v + 35)) - 1); //the probability of a closed gate to open
    bm = 4 * exp(-(v + 60) / 18); //the probability of an open gate to be closed.
    ah = 0.07*exp(-(v + 58) / 20);
    bh = 1 / (exp(-0.1*(v + 28)) + 1);
    an = -0.01*(v + 34) / (exp(-0.1*(v + 34)) - 1);
    bn = 0.125 * exp(-(v + 44) / 80);

    //second stage -> gate states

    dh = phi * ( ah*(1 - h) - bh*h); // inactivation variable h
    dn = phi * ( an*(1 - n) - bn*n); // inward recitifier

    double m_new = am / (am + bm);   //activation variable 
    double h_new = dh*dt + h;
    double n_new = dn*dt + n;

    mvector.push_back(m_new);
    hvector.push_back(h_new); 
    nvector.push_back(n_new);

    // third stage -> currents

    double IL = gL*(v - EL);
    double INa = gNa * pow(m, 3) * h * (v - ENa);
    double IK = gK * pow(n, 4) * (v - EK); // delayed recitifier
    double currents = -INa - IK - IL + Iapp;

    return std::make_tuple(currents);
}


int main() {

             //initialize vectors - later modify with function within function
    double ami = -0.1*(V_init + 35) / (exp(-0.1*(V_init + 35)) - 1);
    double bmi = 4 * exp(-(V_init + 60) / 18);
    double ahi = 0.07*exp(-(V_init + 58) / 20); // 
    double bhi = 1 / (exp(-0.1*(V_init + 28)) + 1);
    double ani = -0.01*(V_init + 34) / (exp(-0.1*(V_init + 34)) - 1);
    double bni = 0.125 * exp(-(V_init + 44) / 80);
    V.push_back(V_init); //V_init
    double m_init = (ami / (ami + bmi));
    double h_init = (ahi / (ahi + bhi)); 
    double n_init = (ani / (ani + bni));

    mvector.push_back(m_init); 
    hvector.push_back(h_init);
    nvector.push_back(n_init);

    for (int i{ 1 }; i <= runtime; ++i) { // 1ms iterations over 3000ms range

        auto ret = gate_probabilities(V[i - 1], mvector[i - 1], hvector[i - 1], nvector[i - 1]);
        double currents;
        currents = std::get<0>(ret);

        double potential = currents*dt + V[i - 1]; //calculate new V.
        V.push_back(potential); //incorporates new V into a vector.
        timevector.push_back( timevector[i - 1] + dt );
    }
}