
        format none;

        function xdot = f (x, t)
            # Zmena hodnoty x za jednotku casu
            xdot = zeros (56,1);
            xdot(1) = 0; # fixed, carb 
xdot(2) = 0; # fixed, icyt 
xdot(3) = 0; # fixed, ocyt 
xdot(4) = 0; # fixed, pps 
xdot(5) = 0; # fixed, liq 
xdot(6) = 0; # fixed, D_C 
xdot(7) = 0; # fixed, D_H 
xdot(8) = 0; # fixed, K_3 
xdot(9) = 0; # fixed, K_4 
xdot(10) = 0; # fixed, K_ba 
xdot(11) = 0; # fixed, K_c 
xdot(12) = 0; # fixed, K_ca 
xdot(13) = 0; # fixed, K_mc 
xdot(14) = 0; # fixed, K_r 
xdot(15) = 0; # fixed, K_t 
xdot(16) = 0; # fixed, N 
xdot(17) = 0; # fixed, P_b 
xdot(18) = 0; # fixed, P_cc 
xdot(19) = 0; # fixed, P_pH 
xdot(20) = 0; # fixed, r_1 
xdot(21) = 0; # fixed, r_2 
xdot(22) = 0; # fixed, r_3 
xdot(23) = 0; # fixed, r_4 
xdot(24) = 0; # fixed, r_5 
xdot(25) = 0; # fixed, r_6 
xdot(26) = 0; # fixed, V_ca 
xdot(27) = 0; # fixed, V_cvolume 
xdot(28) = 0; # fixed, V_e 
xdot(29) = 0; # fixed, V_p 
xdot(30) = 0; # fixed, V_t 
xdot(38) = 0; # fixed, pi 
xdot(39) = 0; # fixed, zero 
xdot(55) = 0; # fixed, HCO3_m_liq 
xdot(56) = 0; # fixed, CO2_liq 
xdot(31) = x(31) = 4*x(38)*x(23)^2; # assignment, A_r_4  FROM INIT EXPRESSION!!! 
xdot(32) = x(32) = 4*x(38)*x(22)^2; # assignment, A_r_3  FROM INIT EXPRESSION!!! 
xdot(33) = x(33) = x(17) / 100; # assignment, P_tH  FROM INIT EXPRESSION!!! 
xdot(34) = x(34) = x(7) / 80; # assignment, D_xC  FROM INIT EXPRESSION!!! 
xdot(35) = x(35) = 0.09*x(34); # assignment, D_xH  FROM INIT EXPRESSION!!! 
xdot(36) = x(36) = 1.18*x(26); # assignment, V_ba  FROM INIT EXPRESSION!!! 
xdot(37) = x(37) = 3.5 * x(29); # assignment, V_cmax  FROM INIT EXPRESSION!!! 
xdot(40) = x(40) = (x(30)*x(53))/(x(15) + x(53)); # assignment, v9  FROM INIT EXPRESSION!!! 
xdot(41) = x(41) = x(8)*x(51)*x(28); # assignment, v5  FROM INIT EXPRESSION!!! 
xdot(42) = x(42) = x(17)*x(31)*(x(54) - x(52)) ; # assignment, v15  FROM INIT EXPRESSION!!! 
xdot(43) = x(43) = (x(29)*x(48))/(x(11) + x(48)) ; # assignment, v8  FROM INIT EXPRESSION!!! 
xdot(44) = x(44) = (x(36)*x(12)*x(47))/(x(10)*x(12)+x(12)*x(47)+x(10)*x(48)) ; # assignment, v1  FROM INIT EXPRESSION!!! 
xdot(45) = x(45) = (x(26)*x(10)*x(48))/(x(12)*x(10)+x(10)*x(48)+x(12)*x(47)) ; # assignment, v2  FROM INIT EXPRESSION!!! 
xdot(46) = x(46) = (6*4*x(38)*x(34)*x(20)*x(21)*(x(48) - x(50)))/(x(21)-x(20)) ; # assignment, v11  FROM INIT EXPRESSION!!! 
xdot(47) =  - ((x(36)*x(12)*x(47))/(x(10)*x(12)+x(12)*x(47)+x(10)*x(48)))  + ((x(16)*4*x(38)*x(35)*x(20)*x(21)*(x(49) - x(47)))/(x(21)-x(20)))  + ((x(26)*x(10)*x(48))/(x(12)*x(10)+x(10)*x(48)+x(12)*x(47))) ; # reaction, HCO3_m_carb 
xdot(48) =  + ((x(36)*x(12)*x(47))/(x(10)*x(12)+x(12)*x(47)+x(10)*x(48)))  - ((x(29)*x(48))/(x(11) + x(48)))  - ((x(26)*x(10)*x(48))/(x(12)*x(10)+x(10)*x(48)+x(12)*x(47)))  - ((x(16)*4*x(38)*x(34)*x(20)*x(21)*(x(48) - x(50)))/(x(21)-x(20))) ; # reaction, CO2_carb 
xdot(49) =  - (x(8)*x(49)*x(27))  - ((x(16)*4*x(38)*x(35)*x(20)*x(21)*(x(49) - x(47)))/(x(21)-x(20)))  + (x(9)*x(50)*x(27))  + (x(33)*x(32)*(x(51) - x(49))  ) ; # reaction, HCO3_m_icyt 
xdot(50) =  + (x(8)*x(49)*x(27))  - (x(9)*x(50)*x(27))  + ((x(16)*4*x(38)*x(34)*x(20)*x(21)*(x(48) - x(50)))/(x(21)-x(20)))  - (x(18)*x(32)*(x(50) - x(52))  ) ; # reaction, CO2_cyt_icyt 
xdot(51) =  + ((x(37)*x(52))/(x(13) + x(52)))  - (x(8)*x(51)*x(28))  - (x(19)*x(31)*(x(51) - x(53))  )  - ((x(30)*x(51))/(x(14) + x(51)))  + ((x(30)*x(53))/(x(15) + x(53)))  + (x(9)*x(52)*x(28))  - (x(33)*x(32)*(x(51) - x(49))  ) ; # reaction, HCO3_m_ocyt 
xdot(52) =  - ((x(37)*x(52))/(x(13) + x(52)))  + (x(8)*x(51)*x(28))  - (x(9)*x(52)*x(28))  + (x(18)*x(32)*(x(50) - x(52))  )  + (x(17)*x(31)*(x(54) - x(52))  ) ; # reaction, CO2_cyt_ocyt 
xdot(53) =  + (x(19)*x(31)*(x(51) - x(53))  )  + ((4*x(38)*x(7)*x(24)*x(25)*(x(55)-x(53)))/(x(25)-x(24)))  + ((x(30)*x(51))/(x(14) + x(51)))  - ((x(30)*x(53))/(x(15) + x(53))) ; # reaction, HCO3_m_pps 
xdot(54) =  + ((4*x(38)*x(6)*x(24)*x(25)*(x(56)-x(54)))/(x(25)-x(24)))  - (x(17)*x(31)*(x(54) - x(52))  ) ; # reaction, CO2_pps 
xdot(31) = 0; # end
xdot(32) = 0; # end
xdot(33) = 0; # end
xdot(34) = 0; # end
xdot(35) = 0; # end
xdot(36) = 0; # end
xdot(37) = 0; # end
xdot(40) = 0; # end
xdot(41) = 0; # end
xdot(42) = 0; # end
xdot(43) = 0; # end
xdot(44) = 0; # end
xdot(45) = 0; # end
xdot(46) = 0; # end

        endfunction

        # time scale
        t = linspace(0, 300, 1000);

        # default values
        x0 = zeros (56,1);
        x0(1) = 1; # fixed compartment, carb
x0(2) = 1; # fixed compartment, icyt
x0(3) = 1; # fixed compartment, ocyt
x0(4) = 1; # fixed compartment, pps
x0(5) = 1; # fixed compartment, liq
x0(6) = 1.88e+3; # fixed parameter, D_C
x0(7) = 1.15e+3 ; # fixed parameter, D_H
x0(8) = 2.63e-3; # fixed parameter, K_3
x0(9) = 37.2e-3; # fixed parameter, K_4
x0(10) = 1800; # fixed parameter, K_ba
x0(11) = 250; # fixed parameter, K_c
x0(12) = 30311; # fixed parameter, K_ca
x0(13) = 0.4; # fixed parameter, K_mc
x0(14) = 200000; # fixed parameter, K_r
x0(15) = 80; # fixed parameter, K_t
x0(16) = 6; # fixed parameter, N
x0(17) = 3.5e+3; # fixed parameter, P_b
x0(18) = 0.01; # fixed parameter, P_cc
x0(19) = 0.3e-2; # fixed parameter, P_pH
x0(20) = 0.01; # fixed parameter, r_1
x0(21) = 0.2; # fixed parameter, r_2
x0(22) = 1.73; # fixed parameter, r_3
x0(23) = 1.77; # fixed parameter, r_4
x0(24) = 1.87; # fixed parameter, r_5
x0(25) = 3.54; # fixed parameter, r_6
x0(26) = 0.0000000008; # fixed parameter, V_ca
x0(27) = 1.0; # fixed parameter, V_cvolume
x0(28) = 1.0; # fixed parameter, V_e
x0(29) = 8e-12; # fixed parameter, V_p
x0(30) = 1.6e-11 ; # fixed parameter, V_t
x0(38) = 3.14159265359; # fixed parameter, pi
x0(39) = 0; # fixed parameter, zero
x0(47) = 3910; # fixed specie, HCO3_m_carb
x0(48) = 259; # fixed specie, CO2_carb
x0(49) = 8870; # fixed specie, HCO3_m_icyt
x0(50) = 3.7; # fixed specie, CO2_cyt_icyt
x0(51) = 8920; # fixed specie, HCO3_m_ocyt
x0(52) = 3.5; # fixed specie, CO2_cyt_ocyt
x0(53) = 142; # fixed specie, HCO3_m_pps
x0(54) = 3.4; # fixed specie, CO2_pps
x0(55) = 143; # fixed specie, HCO3_m_liq
x0(56) = 3.3; # fixed specie, CO2_liq
x0(31) = 4*x0(38)*x0(23)^2; # assignment parameter, A_r_4
x0(32) = 4*x0(38)*x0(22)^2; # assignment parameter, A_r_3
x0(33) = x0(17) / 100; # assignment parameter, P_tH
x0(34) = x0(7) / 80; # assignment parameter, D_xC
x0(35) = 0.09*x0(34); # assignment parameter, D_xH
x0(36) = 1.18*x0(26); # assignment parameter, V_ba
x0(37) = 3.5 * x0(29); # assignment parameter, V_cmax
x0(40) = (x0(30)*x0(53))/(x0(15) + x0(53)); # assignment parameter, v9
x0(41) = x0(8)*x0(51)*x0(28); # assignment parameter, v5
x0(42) = x0(17)*x0(31)*(x0(54) - x0(52)) ; # assignment parameter, v15
x0(43) = (x0(29)*x0(48))/(x0(11) + x0(48)) ; # assignment parameter, v8
x0(44) = (x0(36)*x0(12)*x0(47))/(x0(10)*x0(12)+x0(12)*x0(47)+x0(10)*x0(48)) ; # assignment parameter, v1
x0(45) = (x0(26)*x0(10)*x0(48))/(x0(12)*x0(10)+x0(10)*x0(48)+x0(12)*x0(47)) ; # assignment parameter, v2
x0(46) = (6*4*x0(38)*x0(34)*x0(20)*x0(21)*(x0(48) - x0(50)))/(x0(21)-x0(20)) ; # assignment parameter, v11


        #lsode_options ( 'absolute tolerance', max (50*eps, 0.5e-28) );
        lsode_options ( 'relative tolerance', max (50*eps, 0.5e-28) ); 

        y = lsode ("f", x0, t);
        
        for i = 1:length(t)
            y(i, 31) = 4*y(i, 38)*y(i, 23)^2; # assignment, A_r_4  FROM INIT EXPRESSION!!! 
y(i, 32) = 4*y(i, 38)*y(i, 22)^2; # assignment, A_r_3  FROM INIT EXPRESSION!!! 
y(i, 33) = y(i, 17) / 100; # assignment, P_tH  FROM INIT EXPRESSION!!! 
y(i, 34) = y(i, 7) / 80; # assignment, D_xC  FROM INIT EXPRESSION!!! 
y(i, 35) = 0.09*y(i, 34); # assignment, D_xH  FROM INIT EXPRESSION!!! 
y(i, 36) = 1.18*y(i, 26); # assignment, V_ba  FROM INIT EXPRESSION!!! 
y(i, 37) = 3.5 * y(i, 29); # assignment, V_cmax  FROM INIT EXPRESSION!!! 
y(i, 40) = (y(i, 30)*y(i, 53))/(y(i, 15) + y(i, 53)); # assignment, v9  FROM INIT EXPRESSION!!! 
y(i, 41) = y(i, 8)*y(i, 51)*y(i, 28); # assignment, v5  FROM INIT EXPRESSION!!! 
y(i, 42) = y(i, 17)*y(i, 31)*(y(i, 54) - y(i, 52)) ; # assignment, v15  FROM INIT EXPRESSION!!! 
y(i, 43) = (y(i, 29)*y(i, 48))/(y(i, 11) + y(i, 48)) ; # assignment, v8  FROM INIT EXPRESSION!!! 
y(i, 44) = (y(i, 36)*y(i, 12)*y(i, 47))/(y(i, 10)*y(i, 12)+y(i, 12)*y(i, 47)+y(i, 10)*y(i, 48)) ; # assignment, v1  FROM INIT EXPRESSION!!! 
y(i, 45) = (y(i, 26)*y(i, 10)*y(i, 48))/(y(i, 12)*y(i, 10)+y(i, 10)*y(i, 48)+y(i, 12)*y(i, 47)) ; # assignment, v2  FROM INIT EXPRESSION!!! 
y(i, 46) = (6*4*y(i, 38)*y(i, 34)*y(i, 20)*y(i, 21)*(y(i, 48) - y(i, 50)))/(y(i, 21)-y(i, 20)) ; # assignment, v11  FROM INIT EXPRESSION!!! 

        endfor
        
        [t',y]
        