Решение системы ОДУ методом Рунге-Кутты 4 порядка с граничными условиями

double fnSR(double Cai, double CanSR, double CajSR){
    double Jup = kup * Cai/(Kup + Cai);
    double Jrefill = krefill * (CanSR - CajSR);
    return ki_nSR*Jup - Jrefill;
}
 
double fjSR(double CanSR, double CajSR, double CaSS){
    double Jrefill = krefill * (CanSR - CajSR);
    double Jrel = krel * (CajSR - CaSS);
    return knSR_jSR*Jrefill - Jrel;
}
 
double fSS(double CajSR, double CaSS, double Cai){
    double Jrel = krel * (CajSR - CaSS);
    double Jdiff = kdiff*(CaSS - Cai);
    return kjSR_SS*Jrel - Jdiff;
}
 
double fi(double CaSS, double Cai){
    double Jup = kup * Cai/(Kup + Cai);
    double Jdiff = kdiff*(CaSS - Cai);
    return ki_SS*Jdiff - Jup;
}
 
void Caclocks::Runge_Kutt(double h){
    double k1nSR = fnSR(Caiprev,CanSRprev,CajSRprev);
    double k1jSR = fjSR(CanSRprev,CajSRprev,CaSSprev);
    double k1SS = fSS(CajSRprev,CaSSprev,Caiprev);
    double k1i = fi(CaSSprev,Caiprev);
    double k2nSR = fnSR(Caiprev + h/2 * k1i,CanSRprev + h/2 * k1nSR,CajSRprev + h/2 * k1jSR);
    double k2jSR = fjSR(CanSRprev + h/2 * k1nSR,CajSRprev + h/2 * k1jSR,CaSSprev + h/2 * k1SS);
    double k2SS = fSS(CajSRprev + h/2 * k1jSR,CaSSprev + h/2 * k1SS,Caiprev + h/2 * k1i);
    double k2i = fi(CaSSprev + h/2 * k1SS,Caiprev + h/2 * k1i);
    double k3nSR = fnSR(Caiprev + h/2 * k2i,CanSRprev + h/2 * k2nSR,CajSRprev + h/2 * k2jSR);
    double k3jSR = fjSR(CanSRprev + h/2 * k2nSR,CajSRprev + h/2 * k2jSR,CaSSprev + h/2 * k2SS);
    double k3SS = fSS(CajSRprev + h/2 * k2jSR,CaSSprev + h/2 * k2SS,Caiprev + h/2 * k2i);
    double k3i = fi(CaSSprev + h/2 * k2SS,Caiprev + h/2 * k2i);
    double k4nSR = fnSR(Caiprev + h * k3i,CanSRprev + h * k3nSR,CajSRprev + h * k3jSR);
    double k4jSR = fjSR(CanSRprev + h * k3nSR,CajSRprev + h * k3jSR,CaSSprev + h * k3SS);
    double k4SS = fSS(CajSRprev + h * k3jSR,CaSSprev + h * k3SS,Caiprev + h * k3i);
    double k4i = fi(CaSSprev + h * k3SS,Caiprev + h * k3i);
    CanSRcur = CanSRprev + h/6 * (k1nSR + 2*k2nSR + 2*k3nSR + k4nSR);
    CajSRcur = CajSRprev + h/6 * (k1jSR + 2*k2jSR + 2*k3jSR + k4jSR);
    CaSScur = CaSSprev + h/6 * (k1SS + 2*k2SS + 2*k3SS + k4SS);
    Caicur = Caiprev + h/6 * (k1i + 2*k2i + 2*k3i + k4i);
}