def Battery_Func(t, SV, SV_dot, sw):
    '================================================================='
    '==========================INITIALIZE============================='
    '================================================================='
    print(t)
    nSV = len(SV)
    res = np.zeros([nSV])
    offset_vec = sep.offsets
    "        anode = an.obj['electrode']\n        anode_s = an.obj['surf']\n        elyte = an.obj['elyte']\n        cathode = cat.obj['electrode']\n        cathode_s = cat.obj['surf']"
    nsp_an = anode.n_species
    nsp_cat = cathode.n_species
    F = ct.faraday
    R = ct.gas_constant
    T = Inputs.T
    '================================================================='
    '============================ANODE================================'
    '================================================================='
    offset = an.offsets
    ptr = an.ptr
    j = 0
    N_io_m = 0
    i_io_m = 0
    i_el_m = an.i_ext
    X_an_1 = SV[(offset[j] + ptr['X_ed'][(- 1)])]
    rho_k_elyte_1 = SV[(offset[j] + ptr['rho_k_elyte'])]
    phi_elec_an_1 = SV[(offset[j] + ptr['Phi_ed'])]
    phi_elec_elyte_1 = (phi_elec_an_1 - SV[(offset[j] + ptr['Phi_dl'])])
    anode.X = [X_an_1, (1 - X_an_1)]
    anode.electric_potential = phi_elec_an_1
    conductor.electric_potential = phi_elec_an_1
    elyte.Y = (rho_k_elyte_1 / np.sum(rho_k_elyte_1))
    elyte.electric_potential = phi_elec_elyte_1
    sdot_1 = anode_s.net_production_rates
    j = 1
    offset = int(offset_vec[j])
    X_an_2 = SV[(offset + an.ptr['X_ed'][(- 1)])]
    rho_k_elyte_2 = SV[(offset + an.ptr['rho_k_elyte'])]
    phi_elec_an_2 = SV[(offset + an.ptr['Phi_ed'])]
    phi_elec_elyte_2 = (phi_elec_an_2 - SV[(offset + an.ptr['Phi_dl'])])
    anode.X = [X_an_2, (1 - X_an_2)]
    conductor.electric_potential = phi_elec_an_2
    anode.electric_potential = phi_elec_an_2
    elyte.Y = (rho_k_elyte_2 / np.sum(rho_k_elyte_2))
    elyte.electric_potential = phi_elec_elyte_2
    sdot_2 = anode_s.net_production_rates
    j = 0
    offset = int(offset_vec[j])
    i_el_p = ((an.sigma_eff_ed * (phi_elec_an_1 - phi_elec_an_2)) * an.dyInv)
    N_io_p = ((((- an.u_Li_elyte) * elyte.density_mole) * (((R * T) * (rho_k_elyte_2 - rho_k_elyte_1)) + (F * (phi_elec_elyte_2 - phi_elec_elyte_1)))) * an.dyInv)
    i_io_p = (np.dot(N_io_p, Inputs.z_k_elyte) * F)
    i_Far_1 = (((sdot_1[an.ptr['iFar']] * F) * an.A_surf) / an.dyInv)
    X_Li = (1 - SV[(offset + an.ptr['X_ed'])])
    DiffFlux = np.zeros([(an.nshells + 1)])
    DiffFlux[1:(- 1)] = ((an.D_Li_ed * (X_Li[0:(- 1)] - X_Li[1:])) / an.dr)
    DiffFlux[(- 1)] = (sdot_1[0] / anode.density_mole)
    k_m = (np.arange(0, an.nshells) / an.nshells)
    k_p = (np.arange(1, (an.nshells + 1)) / an.nshells)
    'Calculate the change in X_C6 in the particle interior.\n            Note that the DiffFlux is the diffusion of lithium\n            toward the particle surface, and that diffusion of Li\n            into the shell decreases the amount of C6. The fluxes\n            must be scaled by the shell interface surface area\n            relative to the total particle surface area'
    res[(offset + an.ptr['X_ed'])] = (SV_dot[(offset + an.ptr['X_ed'])] - (((((DiffFlux[1:] * (k_p ** 2)) - (DiffFlux[0:(- 1)] * (k_m ** 2))) * an.A_surf) / an.eps_ed) / an.V_shell))
    'Change in electrolyte_composition'
    res[(offset + an.ptr['rho_k_elyte'])] = (SV_dot[(offset + an.ptr['rho_k_elyte'])] - (((((N_io_m - N_io_p) * an.dyInv) + (sdot_1[nsp_an] * an.A_surf)) / elyte.density_mole) / an.eps_elyte))
    'Double-layer voltage'
    res[(offset + an.ptr['Phi_dl'])] = (SV_dot[(offset + an.ptr['Phi_dl'])] - (((((i_Far_1 + i_el_m) - i_el_p) * an.dyInv) / an.C_dl) / an.A_surf))
    'Algebraic equation for ANODE electric potential boundary condition'
    res[(offset + an.ptr['Phi_ed'])] = SV[(offset + an.ptr['Phi_ed'])]
    '============================ANODE================================'
    'INTERIOR NODES'
    for j in np.arange(2, an.npoints):
        N_io_m = N_io_p
        i_io_m = i_io_p
        i_el_m = i_el_p
        X_an_1 = X_an_2
        rho_k_elyte_1 = rho_k_elyte_2
        phi_elec_an_1 = phi_elec_an_2
        phi_elec_elyte_1 = phi_elec_elyte_2
        sdot_1 = sdot_2
        offset = int(an.offsets[j])
        X_an_2 = SV[(offset + an.ptr['X_ed'][(- 1)])]
        rho_k_elyte_2 = SV[(offset + an.ptr['rho_k_elyte'])]
        phi_elec_an_2 = SV[(offset + an.ptr['Phi_ed'])]
        phi_elec_elyte_2 = (phi_elec_an_2 - SV[(offset + an.ptr['Phi_dl'])])
        anode.X = [X_an_2, (1 - X_an_2)]
        anode.electric_potential = phi_elec_an_2
        conductor.electric_potential = phi_elec_an_2
        elyte.Y = (rho_k_elyte_2 / np.sum(rho_k_elyte_2))
        elyte.electric_potential = phi_elec_elyte_2
        sdot_2 = anode_s.net_production_rates
        offset = int(an.offsets[(j - 1)])
        i_el_p = ((an.sigma_eff_ed * (phi_elec_an_1 - phi_elec_an_2)) * an.dyInv)
        N_io_p = ((((- an.u_Li_elyte) * elyte.density_mole) * (((R * T) * (rho_k_elyte_2 - rho_k_elyte_1)) + (F * (phi_elec_elyte_2 - phi_elec_elyte_1)))) * an.dyInv)
        i_io_p = (np.dot(N_io_p, Inputs.z_k_elyte) * F)
        i_Far_1 = (((sdot_1[an.ptr['iFar']] * F) * an.A_surf) / an.dyInv)
        X_Li = (1 - SV[(offset + an.ptr['X_ed'])])
        DiffFlux = np.zeros([(an.nshells + 1)])
        DiffFlux[1:(- 1)] = ((an.D_Li_ed * (X_Li[0:(- 1)] - X_Li[1:])) / an.dr)
        DiffFlux[(- 1)] = (sdot_1[0] / anode.density_mole)
        'Calculate the change in X_C6 in the particle interior.'
        res[(offset + an.ptr['X_ed'])] = (SV_dot[(offset + an.ptr['X_ed'])] - (((((DiffFlux[1:] * (k_p ** 2)) - (DiffFlux[0:(- 1)] * (k_m ** 2))) * an.A_surf) / an.eps_ed) / an.V_shell))
        'Change in electrolyte_composition'
        res[(offset + an.ptr['rho_k_elyte'])] = (SV_dot[(offset + an.ptr['rho_k_elyte'])] - (((((N_io_m - N_io_p) * an.dyInv) + (sdot_1[nsp_an] * an.A_surf)) / elyte.density_mole) / an.eps_elyte))
        'Double-layer voltage'
        res[(offset + an.ptr['Phi_dl'])] = (SV_dot[(offset + an.ptr['Phi_dl'])] - (((((i_Far_1 + i_el_m) - i_el_p) * an.dyInv) / an.C_dl) / an.A_surf))
        'Algebraic equation for ANODE electric potential boundary condition'
        res[(offset + an.ptr['Phi_ed'])] = (((i_el_m - i_el_p) + i_io_m) - i_io_p)
    '============================ANODE================================'
    'Separator boundary'
    N_io_m = N_io_p
    i_io_m = i_io_p
    i_el_m = i_el_p
    X_an_1 = X_an_2
    rho_k_elyte_1 = rho_k_elyte_2
    phi_elec_an_1 = phi_elec_an_2
    phi_elec_elyte_1 = phi_elec_elyte_2
    sdot_1 = sdot_2
    i_el_p = 0
    j = (an.npoints - 1)
    offset = int(an.offsets[j])
    i_Far_1 = (((sdot_1[an.ptr['iFar']] * F) * an.A_surf) / an.dyInv)
    i_io_p = an.i_ext
    N_io_p = np.zeros_like(N_io_p)
    N_io_p[2] = (i_io_p / F)
    X_Li = (1 - SV[(offset + an.ptr['X_ed'])])
    DiffFlux = np.zeros([(an.nshells + 1)])
    DiffFlux[1:(- 1)] = ((an.D_Li_ed * (X_Li[0:(- 1)] - X_Li[1:])) / an.dr)
    DiffFlux[(- 1)] = (sdot_1[0] / anode.density_mole)
    'Calculate the change in X_C6 in the particle interior.'
    res[(offset + an.ptr['X_ed'])] = (SV_dot[(offset + an.ptr['X_ed'])] - (((((DiffFlux[1:] * (k_p ** 2)) - (DiffFlux[0:(- 1)] * (k_m ** 2))) * an.A_surf) / an.eps_ed) / an.V_shell))
    'Change in electrolyte_composition'
    res[(offset + an.ptr['rho_k_elyte'])] = (SV_dot[(offset + an.ptr['rho_k_elyte'])] - (((((N_io_m - N_io_p) * an.dyInv) + (sdot_1[nsp_an] * an.A_surf)) / elyte.density_mole) / an.eps_elyte))
    'Double-layer voltage'
    res[(offset + an.ptr['Phi_dl'])] = (SV_dot[(offset + an.ptr['Phi_dl'])] - (((((i_Far_1 + i_el_m) - i_el_p) * an.dyInv) / an.C_dl) / an.A_surf))
    'Algebraic equation for ANODE electric potential boundary condition'
    res[(offset + an.ptr['Phi_ed'])] = (((i_el_m - i_el_p) + i_io_m) - i_io_p)
    '================================================================='
    '==========================SEPARATOR=============================='
    '================================================================='
    '================================================================='
    '===========================CATHODE==============================='
    '================================================================='
    '=========================CATHODE============================='
    'current collector boundary'
    return res
=================================================================
def Battery_Func(t, SV, SV_dot, sw):
    
    '==========================INITIALIZE============================='
    
    print(t)
    nSV = len(SV)
    res = np.zeros([nSV])
    offset_vec = sep.offsets
    "        anode = an.obj['electrode']\n        anode_s = an.obj['surf']\n        elyte = an.obj['elyte']\n        cathode = cat.obj['electrode']\n        cathode_s = cat.obj['surf']"
    nsp_an = anode.n_species
    nsp_cat = cathode.n_species
    F = ct.faraday
    R = ct.gas_constant
    T = Inputs.T
    
    '============================ANODE================================'
    
    offset = an.offsets
    ptr = an.ptr
    j = 0
    N_io_m = 0
    i_io_m = 0
    i_el_m = an.i_ext
    X_an_1 = SV[(offset[j] + ptr['X_ed'][(- 1)])]
    rho_k_elyte_1 = SV[(offset[j] + ptr['rho_k_elyte'])]
    phi_elec_an_1 = SV[(offset[j] + ptr['Phi_ed'])]
    phi_elec_elyte_1 = (phi_elec_an_1 - SV[(offset[j] + ptr['Phi_dl'])])
    anode.X = [X_an_1, (1 - X_an_1)]
    anode.electric_potential = phi_elec_an_1
    conductor.electric_potential = phi_elec_an_1
    elyte.Y = (rho_k_elyte_1 / np.sum(rho_k_elyte_1))
    elyte.electric_potential = phi_elec_elyte_1
    sdot_1 = anode_s.net_production_rates
    j = 1
    offset = int(offset_vec[j])
    X_an_2 = SV[(offset + an.ptr['X_ed'][(- 1)])]
    rho_k_elyte_2 = SV[(offset + an.ptr['rho_k_elyte'])]
    phi_elec_an_2 = SV[(offset + an.ptr['Phi_ed'])]
    phi_elec_elyte_2 = (phi_elec_an_2 - SV[(offset + an.ptr['Phi_dl'])])
    anode.X = [X_an_2, (1 - X_an_2)]
    conductor.electric_potential = phi_elec_an_2
    anode.electric_potential = phi_elec_an_2
    elyte.Y = (rho_k_elyte_2 / np.sum(rho_k_elyte_2))
    elyte.electric_potential = phi_elec_elyte_2
    sdot_2 = anode_s.net_production_rates
    j = 0
    offset = int(offset_vec[j])
    i_el_p = ((an.sigma_eff_ed * (phi_elec_an_1 - phi_elec_an_2)) * an.dyInv)
    N_io_p = ((((- an.u_Li_elyte) * elyte.density_mole) * (((R * T) * (rho_k_elyte_2 - rho_k_elyte_1)) + (F * (phi_elec_elyte_2 - phi_elec_elyte_1)))) * an.dyInv)
    i_io_p = (np.dot(N_io_p, Inputs.z_k_elyte) * F)
    i_Far_1 = (((sdot_1[an.ptr['iFar']] * F) * an.A_surf) / an.dyInv)
    X_Li = (1 - SV[(offset + an.ptr['X_ed'])])
    DiffFlux = np.zeros([(an.nshells + 1)])
    DiffFlux[1:(- 1)] = ((an.D_Li_ed * (X_Li[0:(- 1)] - X_Li[1:])) / an.dr)
    DiffFlux[(- 1)] = (sdot_1[0] / anode.density_mole)
    k_m = (np.arange(0, an.nshells) / an.nshells)
    k_p = (np.arange(1, (an.nshells + 1)) / an.nshells)
    'Calculate the change in X_C6 in the particle interior.\n            Note that the DiffFlux is the diffusion of lithium\n            toward the particle surface, and that diffusion of Li\n            into the shell decreases the amount of C6. The fluxes\n            must be scaled by the shell interface surface area\n            relative to the total particle surface area'
    res[(offset + an.ptr['X_ed'])] = (SV_dot[(offset + an.ptr['X_ed'])] - (((((DiffFlux[1:] * (k_p ** 2)) - (DiffFlux[0:(- 1)] * (k_m ** 2))) * an.A_surf) / an.eps_ed) / an.V_shell))
    'Change in electrolyte_composition'
    res[(offset + an.ptr['rho_k_elyte'])] = (SV_dot[(offset + an.ptr['rho_k_elyte'])] - (((((N_io_m - N_io_p) * an.dyInv) + (sdot_1[nsp_an] * an.A_surf)) / elyte.density_mole) / an.eps_elyte))
    'Double-layer voltage'
    res[(offset + an.ptr['Phi_dl'])] = (SV_dot[(offset + an.ptr['Phi_dl'])] - (((((i_Far_1 + i_el_m) - i_el_p) * an.dyInv) / an.C_dl) / an.A_surf))
    'Algebraic equation for ANODE electric potential boundary condition'
    res[(offset + an.ptr['Phi_ed'])] = SV[(offset + an.ptr['Phi_ed'])]
    '============================ANODE================================'
    'INTERIOR NODES'
    for j in np.arange(2, an.npoints):
        N_io_m = N_io_p
        i_io_m = i_io_p
        i_el_m = i_el_p
        X_an_1 = X_an_2
        rho_k_elyte_1 = rho_k_elyte_2
        phi_elec_an_1 = phi_elec_an_2
        phi_elec_elyte_1 = phi_elec_elyte_2
        sdot_1 = sdot_2
        offset = int(an.offsets[j])
        X_an_2 = SV[(offset + an.ptr['X_ed'][(- 1)])]
        rho_k_elyte_2 = SV[(offset + an.ptr['rho_k_elyte'])]
        phi_elec_an_2 = SV[(offset + an.ptr['Phi_ed'])]
        phi_elec_elyte_2 = (phi_elec_an_2 - SV[(offset + an.ptr['Phi_dl'])])
        anode.X = [X_an_2, (1 - X_an_2)]
        anode.electric_potential = phi_elec_an_2
        conductor.electric_potential = phi_elec_an_2
        elyte.Y = (rho_k_elyte_2 / np.sum(rho_k_elyte_2))
        elyte.electric_potential = phi_elec_elyte_2
        sdot_2 = anode_s.net_production_rates
        offset = int(an.offsets[(j - 1)])
        i_el_p = ((an.sigma_eff_ed * (phi_elec_an_1 - phi_elec_an_2)) * an.dyInv)
        N_io_p = ((((- an.u_Li_elyte) * elyte.density_mole) * (((R * T) * (rho_k_elyte_2 - rho_k_elyte_1)) + (F * (phi_elec_elyte_2 - phi_elec_elyte_1)))) * an.dyInv)
        i_io_p = (np.dot(N_io_p, Inputs.z_k_elyte) * F)
        i_Far_1 = (((sdot_1[an.ptr['iFar']] * F) * an.A_surf) / an.dyInv)
        X_Li = (1 - SV[(offset + an.ptr['X_ed'])])
        DiffFlux = np.zeros([(an.nshells + 1)])
        DiffFlux[1:(- 1)] = ((an.D_Li_ed * (X_Li[0:(- 1)] - X_Li[1:])) / an.dr)
        DiffFlux[(- 1)] = (sdot_1[0] / anode.density_mole)
        'Calculate the change in X_C6 in the particle interior.'
        res[(offset + an.ptr['X_ed'])] = (SV_dot[(offset + an.ptr['X_ed'])] - (((((DiffFlux[1:] * (k_p ** 2)) - (DiffFlux[0:(- 1)] * (k_m ** 2))) * an.A_surf) / an.eps_ed) / an.V_shell))
        'Change in electrolyte_composition'
        res[(offset + an.ptr['rho_k_elyte'])] = (SV_dot[(offset + an.ptr['rho_k_elyte'])] - (((((N_io_m - N_io_p) * an.dyInv) + (sdot_1[nsp_an] * an.A_surf)) / elyte.density_mole) / an.eps_elyte))
        'Double-layer voltage'
        res[(offset + an.ptr['Phi_dl'])] = (SV_dot[(offset + an.ptr['Phi_dl'])] - (((((i_Far_1 + i_el_m) - i_el_p) * an.dyInv) / an.C_dl) / an.A_surf))
        'Algebraic equation for ANODE electric potential boundary condition'
        res[(offset + an.ptr['Phi_ed'])] = (((i_el_m - i_el_p) + i_io_m) - i_io_p)
    '============================ANODE================================'
    'Separator boundary'
    N_io_m = N_io_p
    i_io_m = i_io_p
    i_el_m = i_el_p
    X_an_1 = X_an_2
    rho_k_elyte_1 = rho_k_elyte_2
    phi_elec_an_1 = phi_elec_an_2
    phi_elec_elyte_1 = phi_elec_elyte_2
    sdot_1 = sdot_2
    i_el_p = 0
    j = (an.npoints - 1)
    offset = int(an.offsets[j])
    i_Far_1 = (((sdot_1[an.ptr['iFar']] * F) * an.A_surf) / an.dyInv)
    i_io_p = an.i_ext
    N_io_p = np.zeros_like(N_io_p)
    N_io_p[2] = (i_io_p / F)
    X_Li = (1 - SV[(offset + an.ptr['X_ed'])])
    DiffFlux = np.zeros([(an.nshells + 1)])
    DiffFlux[1:(- 1)] = ((an.D_Li_ed * (X_Li[0:(- 1)] - X_Li[1:])) / an.dr)
    DiffFlux[(- 1)] = (sdot_1[0] / anode.density_mole)
    'Calculate the change in X_C6 in the particle interior.'
    res[(offset + an.ptr['X_ed'])] = (SV_dot[(offset + an.ptr['X_ed'])] - (((((DiffFlux[1:] * (k_p ** 2)) - (DiffFlux[0:(- 1)] * (k_m ** 2))) * an.A_surf) / an.eps_ed) / an.V_shell))
    'Change in electrolyte_composition'
    res[(offset + an.ptr['rho_k_elyte'])] = (SV_dot[(offset + an.ptr['rho_k_elyte'])] - (((((N_io_m - N_io_p) * an.dyInv) + (sdot_1[nsp_an] * an.A_surf)) / elyte.density_mole) / an.eps_elyte))
    'Double-layer voltage'
    res[(offset + an.ptr['Phi_dl'])] = (SV_dot[(offset + an.ptr['Phi_dl'])] - (((((i_Far_1 + i_el_m) - i_el_p) * an.dyInv) / an.C_dl) / an.A_surf))
    'Algebraic equation for ANODE electric potential boundary condition'
    res[(offset + an.ptr['Phi_ed'])] = (((i_el_m - i_el_p) + i_io_m) - i_io_p)
    
    '==========================SEPARATOR=============================='
    
    
    '===========================CATHODE==============================='
    
    '=========================CATHODE============================='
    'current collector boundary'
    return res