diff --git a/Model.py b/Model.py
new file mode 100644
index 0000000..627e59c
--- /dev/null
+++ b/Model.py
@@ -0,0 +1,407 @@
+# -*- coding: utf-8 -*-
+
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.integrate import ode
+
+# %% Constants
+
+R = 8.314  # ideal gas constant (J/(mol*K))
+F = 96.485  # coulomb_per_millimole, Faraday constant
+T = 298  # room temperature (K)
+RT = R * T
+
+
+# %% model
+class Model:
+    def __init__(self):
+
+        self.A_cap = 1.534e-4  # cm^2, Capacitive membrane area
+        self.C_mem = 1.0  # uF/cm2 Specific membrane capacitance
+        self.V_myo = 25.84e-6  # uL, Myoplasmic volume
+        self.current2flux = self.A_cap * self.C_mem / 2 / F  #
+        self.period = 1000  # ms, pulse period
+        self.V_mem_rest = -80.0  # mV, resting membrane potential
+        self.Nai = 11000  # uM, Myoplasmic Na+ concentration
+        self.Nao = 150000  # uM, Extracellular Na+ concentration
+        self.eta = 0.35  # Controls voltage dependance of Na/Ca2+ exchange
+        self.km_Na = 87500  # uM, Na+ half-saturation constant for Na+/Ca2+ exchange
+        self.k_sat = (
+            0.1  # Na+/Ca2+ exchange saturation factor at very negative potentials
+        )
+        self.km_Ca = (1380,)  # uM, Ca2+ half-saturation constant for Na+/Ca2+ exchange
+        self.k_Na_Ca = 292.8  # pA/pF, Scaling factor of Na2+/Ca2+ exchange
+        self.I_max_pCa = 1.0  # pA/pF, Maximum Ca2+ pump current
+        # self.Cai = 0.1  # uM, Cytoplasmic Ca2+ concentration fixed at this point
+        self.LTRPN_tot = (
+            70.0  # uM, Total myoplasmic troponin low-affinity site concentration
+        )
+        self.HTRPN_tot = (
+            140.0  # uM, Total myoplasmic troponin high-affinity site concentration
+        )
+        # self.LTRPNCa = 11.2684 # uM, Concentration Ca2+ bound low-affinity troponin-binding sites
+        # self.HTRPNCa = 125.290 # uM, Concentration Ca2+ bound high-affinity troponin-binding sites
+        self.k_htrpn_positive = (
+            0.00237  # uM^(-1)/ms, Ca2+ on rate const. for troponin high-affinity sites
+        )
+        self.k_htrpn_negative = (
+            3.2e-5  # ms, Ca2+ off rate const. for troponin high-affinity sites
+        )
+        self.k_ltrpn_positive = (
+            0.0327  # uM^(-1)/ms, Ca2+ on rate const. for troponin low-affinity sites
+        )
+        self.k_ltrpn_negative = (
+            0.0196  # ms, Ca2+ off rate const. for troponin low-affinity sites
+        )
+        self.CMDN_tot = 50  # uM, Total myoplasmic calmodulin concentration
+        self.Km_CMDN = 0.238  # uM, Ca2 half-saturation constant for calmodulin
+        self.nu_1 = 4.5  # 1/ms, Maximum RyR channel Ca2+ permeability Ca2+ leak rate constant from the NSR
+        self.nu_2 = 1.74e-5  # ms^(-1), Ca2+ leak rate const. from the NSR
+        self.nu_3 = 0.45  # uM/ms, SR Ca2+ -ATPase maximum pupmp rate
+        self.Km_up = 0.5  # uM, Half-saturation constant for SR Ca2+ -ATPase pump
+        self.km_p_ca = 0.5  # uM, Ca2+ half-saturation constant for Ca2+ pump current
+        # self.Ca_NSR = 1299.50 # uM,NSR Ca2+ concentration
+        self.tau_xfer = 8.0  # ms, Time constant for transfer from subspace to myoplasm
+        self.K_pc_max = (
+            0.23324  # 1/ms, Maximum time constant for Ca2+-induced inactivation
+        )
+        self.K_pc_half = (
+            20.0  # uM, Half-saturation constant for Ca2+-induced inactivation
+        )
+        self.Kpcb = 0.0005  # 1/ms, Voltage-insensitive rate constant for inactivation
+        self.km_Ca = 1380  # uM, Ca2+ half-saturation constant for Na+/Ca2+ exchange
+        self.Gcab = 0.000367  # mS/uF
+        self.Cao = 1130.0  # uM, Ca2+ outside the cell
+        self.gCaL = (
+            0.1729  # mS/uF, Specific maximum conductivity for L-type Ca2+ channel
+        )
+        self.ECaL = 43.0  # mV, Reversal potential for L-type Ca2⫹ channel, kas arvutame voi jaab const?
+        self.V_ss = 1.485e-9  # uL, Dyadic aka subspace volume
+        # self.Ca_JSR = 1299.50 # uM, JSR Ca2+ concentration
+        self.F_tot = 25  # uM, total concentration of Fluo-4
+        self.k_on = 0.1  # 1/uM * 1/ms, Fluo-4 reaction rate constant
+        self.k_off = 0.11  # 1/ms, Fluo-4 reaction rate constant 2
+        self.V_NSR = 2.098e-6  # ul, Network SR volume
+        self.tau_tr = 20  # ms, Time const for transfer from NSR to JSR
+        self.CSQN_tot = 15000.0  # uM, total junctional SR calsequestrin concentration
+        self.Km_CSQN = 800.0  # uM, Ca2 half-saturation constant for calsequestrin
+        self.k_a_positive = 0.006075  # uM^(-4)/ms, RyR Pc1 - Po1 rate constant
+        self.k_a_negative = 0.07125  # 1/ms, RyR Po1 - Pc1 rate constant
+        self.k_b_positive = 0.00405  # uM^(-3)/ms, RyR Po1 - Po2 rate constant
+        self.k_b_negative = 0.965  # 1/ms, RyR Po2 - Po1 rate constant
+        self.k_c_positive = 0.009  # 1/ms, RyR Po1 - Pc2 rate constant
+        self.k_c_negative = 0.0008  # 1/ms, RyR Pc2 - Po1 rate constant
+        self.n_ryr = 4  # RyR Ca2+ cooperativity parameter Pc1 - Po1
+        self.m_ryr = 3  # RyR Ca2+ cooperativity parameter Po1 - Po2
+        self.I_CaL_max = 7.0  # pA/pF, normalization constant for L-type Ca2+ current
+        self.V_JSR = 0.12e-6  # ul, Junctional SR volume
+
+    def ode_system(self, t, states):
+
+        (
+            Cass,
+            Os,
+            C2,
+            C3,
+            C4,
+            I1,
+            I2,
+            I3,
+            LTRPNCa,
+            HTRPNCa,
+            Ca_NSR,
+            Ca_JSR,
+            P_RyR,
+            P_O1,
+            P_O2,
+            P_C2,
+            Cai,
+            FCa,
+        ) = states
+
+        V = self.mem_potential(t)
+        VF = V * F
+
+        Kpcf = 13.0 * (1 - np.exp(-((V + 14.5) ** 2) / 100))
+
+        alpha = (  # dimensionless parameter for L-type Ca2+ channel closed states
+            0.4
+            * np.exp((V + 12.0) / 10)
+            * (
+                1
+                + 0.7 * np.exp((-((V + 40.0) ** 2)) / 10.0)
+                - 0.75 * np.exp(-((V + 20.0) ** 2) / 400.0)
+            )
+        ) / (1 + 0.12 * np.exp((V + 12.0) / 10.0))
+
+        beta = 0.05 * np.exp(
+            -(V + 12.0) / 13.0
+        )  # parameter for L-type Ca2+ channel closed states
+
+        gamma = self.K_pc_max * Cass / (self.K_pc_half + Cass)
+
+        C1 = 1 - (Os + C2 + C3 + C4 + I1 + I2 + I3)
+
+        dC2dt = 4 * alpha * C1 - beta * C2 + 2 * beta * C3 - 3 * alpha * C2
+
+        dC3dt = 3 * alpha * C2 - 2 * beta * C3 + 3 * beta * C4 - 2 * alpha * C3
+
+        dC4dt = (
+            2 * alpha * C3
+            - 3 * beta * C4
+            + 4 * beta * Os
+            - alpha * C4
+            + 0.01 * (4 * self.Kpcb * beta * I1 - alpha * gamma * C4)
+            + 0.002 * (4 * beta * I2 - Kpcf * C4)
+            + 4 * beta * self.Kpcb * I3
+            - gamma * Kpcf * C4
+        )
+
+        dI1dt = (
+            gamma * Os
+            - self.Kpcb * I1
+            + 0.001 * (alpha * I3 - Kpcf * I1)
+            + 0.01 * (alpha * gamma * C4 - 4 * beta * self.Kpcb * I1)
+        )
+
+        dI2dt = 0.001 * (Kpcf * Os - alpha * I2)
+        +self.Kpcb * I3 - gamma * I2 + 0.002 * (Kpcf * C4 - 4 * beta * I2)
+
+        dI3dt = (
+            0.001 * (Kpcf * I1 - alpha * I3)
+            + gamma * I2
+            - self.Kpcb * I3
+            + gamma * Kpcf * C4
+            - 4 * beta * self.Kpcb * I3
+        )
+
+        # Bi = (1 + (self.CMDN_tot * self.Km_CMDN)/(self.Km_CMDN+(Cai))**2)**(-1)
+        # Bss = 1 / ((self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cass) ** 2)
+
+        Bi = 1 / (
+            1 + (self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cai) ** 2
+        )  # A6 !!!!!!! korras V8
+
+        Bss = 1 / (
+            1 + (self.CMDN_tot * self.Km_CMDN) / (self.Km_CMDN + Cass) ** 2
+        )  # A7 !!!!!!!!! korras V8
+
+        J_xfer = (Cass - Cai) / self.tau_xfer
+
+        dLTRPNCadt = self.k_ltrpn_positive * Cai * (self.LTRPN_tot - LTRPNCa)
+        -self.k_ltrpn_negative * (LTRPNCa)
+
+        dHTRPNCadt = self.k_htrpn_positive * Cai * (self.HTRPN_tot - HTRPNCa)
+        -self.k_htrpn_negative * (HTRPNCa)
+
+        J_trpn = (
+            self.k_htrpn_positive * Cai * (self.HTRPN_tot - HTRPNCa)
+            - self.k_htrpn_negative * HTRPNCa
+            + self.k_ltrpn_positive
+            * Cai
+            * (self.LTRPN_tot - LTRPNCa)  # !!!!!! korras v8
+            - self.k_ltrpn_negative * LTRPNCa
+        )
+
+        J_up = self.nu_3 * (Cai**2 / (self.Km_up**2 - Cai**2))
+
+        J_tr = (Ca_NSR - Ca_JSR) / self.tau_tr
+
+        J_leak = self.nu_2 * (Ca_NSR - Cai)  # puudu CaNSR diff vorrand
+
+        dCa_NSRdt = (J_up - J_leak) * (self.V_myo / self.V_NSR)
+        -J_tr * (self.V_JSR / self.V_NSR)
+
+        B_JSR = 1 / (1 + (self.CSQN_tot * self.Km_CSQN) / (self.Km_CSQN + Ca_JSR) ** 2)
+
+        J_rel = 0  # self.nu_1*(P_O1+P_O2)*(Ca_JSR-Cass)*P_RyR
+
+        dCa_JSRdt = B_JSR * (J_tr - J_rel)
+
+        b = 1 / (
+            (self.km_Na**3 + self.Nao**3)
+            * (self.km_Ca + self.Cao)
+            * (1 + self.k_sat * np.exp((self.eta - 1) * VF / RT))
+        )
+        p = (
+            np.exp(self.eta * VF / RT) * self.Nai**3 * self.Cao
+            - np.exp((self.eta - 1) * VF / RT) * self.Nao**3 * Cai
+        )
+        I_NaCa = self.k_Na_Ca * b * p
+
+        I_pCa = self.I_max_pCa * ((Cai) ** 2 / ((self.km_p_ca) ** 2 + (Cai) ** 2))
+
+        I_CaL = self.gCaL * Os * (V - self.ECaL)
+
+        dP_RyRdt = -0.04 * P_RyR - 0.1 * (I_CaL / self.I_CaL_max) * np.exp(
+            -((V - 5.0) ** 2) / 648
+        )
+
+        P_C1 = 1 - (P_C2 + P_O1 + P_O2)
+
+        dP_O1dt = self.k_a_positive * (Cass) ** self.n_ryr * P_C1
+        -self.k_a_negative * P_O1 - self.k_a_positive * (Cass) ** self.m_ryr * P_O1
+        +self.k_b_negative * P_O2 - self.k_c_positive * P_O1 + self.k_c_negative * P_C2
+
+        dP_O2dt = (
+            self.k_b_positive * (Cass) ** self.m_ryr * P_O1 - self.k_b_negative * P_O2
+        )
+
+        dP_C2dt = self.k_c_positive * P_O1 - self.k_c_negative * P_C2
+
+        dCassdt = Bss * (
+            -J_xfer * self.V_myo / self.V_ss - I_CaL * self.current2flux / self.V_ss
+        )
+
+        dOsdt = alpha * C4 - 4 * beta * Os + self.Kpcb * I1
+        -gamma * Os + 0.001 * (alpha * I2 - Kpcf * Os)
+
+        ECaN = (R * T) / F * np.log(self.Cao / Cai)
+
+        I_Cab = self.Gcab * (V - ECaN)
+
+        """
+        dCaidt = Bi * (J_leak + J_xfer - J_up - J_trpn
+        - (I_Cab - 2 * I_NaCa + I_pCa) * ((self.A_cap * self.C_mem)/(2 * self.V_myo * F)))
+        """
+
+        # J_rel_caf = self.nu_1 * (self.Ca_JSR - Cass)
+
+        JFCa = self.k_on * (self.F_tot - FCa) * Cai - self.k_off * FCa
+
+        dFCadt = JFCa
+
+        dCaidt = Bi * (
+            J_leak
+            + J_xfer
+            - J_up
+            - J_trpn
+            - JFCa
+            - (I_Cab - 2 * I_NaCa + I_pCa)
+            * ((self.A_cap * self.C_mem) / (2 * self.V_myo * F))
+        )
+
+        # Cass, Os, C2, C3, C4, I1, I2, I3, LTRPNCa, HTRPNCa, Ca_NSR, Ca_JSR, P_RyR, P_O1, P_O2, P_C2, Cai, FCa = states
+
+        return [
+            dCassdt,
+            dOsdt,
+            dC2dt,
+            dC3dt,
+            dC4dt,
+            dI1dt,
+            dI2dt,
+            dI3dt,
+            dLTRPNCadt,
+            dHTRPNCadt,
+            dCa_NSRdt,
+            dCa_JSRdt,
+            dP_RyRdt,
+            dP_O1dt,
+            dP_O2dt,
+            dP_C2dt,
+            dCaidt,
+            dFCadt,
+        ]
+
+    def mem_potential(self, t):
+        v1 = 0
+        t0 = 100
+        t1 = 350
+        n = (t // self.period) * self.period
+        return v1 if t0 + n <= t < t1 + n else self.V_mem_rest
+
+    def solve(self, initial_values, tspan, dt):
+        times = np.arange(*tspan, dt)
+
+        r = ode(self.ode_system)
+        r.set_integrator(
+            "lsoda", method="bdf", atol=1e-07, rtol=1e-07, max_step=0.1, nsteps=500
+        )
+        r.set_initial_value(initial_values, times[0])
+
+        states = np.array([[0.0] * times.size] * len(initial_values))
+        states[:, 0] = initial_values
+
+        for i, t in enumerate(times[1:]):
+            if r.successful():
+                r.integrate(t)
+                states[:, i + 1] = r.y
+            else:
+                break
+
+        V = np.array([self.mem_potential(t) for t in times])
+
+        # plt.plot(times, V)
+        fig = plt.figure()
+        ax1 = fig.add_subplot(141)
+        ax2 = fig.add_subplot(142)
+        ax3 = fig.add_subplot(143)
+        ax4 = fig.add_subplot(144)
+
+        ax1.plot(times, states[0, :], label="Cass")
+        Cai = states[-2, :]
+        ax1.plot(times, (states[0, :] - Cai) / self.tau_xfer, label="Jxfer")
+        ax1.plot(
+            times,
+            -self.gCaL
+            * states[1, :]
+            * (V - self.ECaL)
+            * self.current2flux
+            / self.V_ss
+            / 1000,
+            label="JCaL mM/s",
+        )
+        ax1.legend(frameon=False)
+        ax1.set_xlabel("time [ms]")
+        ax1.set_ylabel(r"$\mu mol/l$")
+
+        ax2.plot(times, self.gCaL * states[1, :] * (V - self.ECaL), label="ICaL ")
+        ax2.legend(frameon=False)
+        ax2.set_xlabel("time [ms]")
+        ax2.set_ylabel("pA/pF")
+
+        ax3.plot(times, Cai, label="Cai ")
+        ax3.legend(frameon=False)
+        ax3.set_xlabel("time [ms]")
+        ax3.set_ylabel(r"$\mu mol/3l$")
+
+        ax4.plot(times, states[-1, :], label="Fca")
+        ax4.legend(frameon=False)
+        ax4.set_xlabel("time [ms]")
+        ax4.set_ylabel(r"$\mu mol/3l$")
+
+        plt.show()
+
+
+if __name__ == "__main__":
+
+    # Cass, Os, C2, C3, C4, I1, I2, I3, LTRPNCa, HTRPNCa, Cai, FCa= states
+    model = Model()
+
+    Cai_0 = 0.11712  # Cai
+    FCa_0 = (model.k_on * model.F_tot * Cai_0) / (model.k_on * Cai_0 + model.k_off)
+    print(FCa_0)
+
+    # Cass, Os, C2, C3, C4, I1, I2, I3, LTRPNCa, HTRPNCa, Ca_NSR, Ca_JSR, P_RyR, P_O1, P_O2, P_C2, Cai, FCa = states
+    initial_values = [
+        Cai_0,  # Cass
+        0.930308e-18,  # Os,
+        0.124216e-3,  # C2,
+        0.578679e-8,  # C3,
+        0.119816e-12,  # C4,
+        0.497923e-18,  # I1,
+        0.345847e-13,  # I2,
+        0.185106e-13,  # I3,
+        11.2684,  # LTRPNCa,
+        125.290,  # HTRPNCa
+        1299.50,  # Ca_NSR
+        1299.50,  # Ca_JSR
+        0.0,  # P_RyR
+        0.149102e-4,  # P_O1
+        0.951726e-10,  # P_02
+        0.167740e-3,  # P_C2
+        Cai_0,  # Cai
+        FCa_0,  # FCa
+    ]
+    model.solve(initial_values=initial_values, tspan=[0, 1000], dt=1.0)