diff --git a/julia/tug/Core/BTCS.jl b/julia/tug/Core/BTCS.jl
index f9d3caf..648aa2a 100644
--- a/julia/tug/Core/BTCS.jl
+++ b/julia/tug/Core/BTCS.jl
@@ -135,38 +135,28 @@ function writeSolutionVector!(sv::Vector{T}, concentrations::Matrix{T}, alphaX::
 end
 
 # BTCS solution for 1D grid
-function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, timestep::T) where {T}
-    length = getCols(grid)
+function BTCS_1D(grid::Grid{T}, bc::Boundary{T}, alpha_left::Matrix{T}, alpha_right::Matrix{T}, timestep::T) where {T}
     sx = timestep / (grid.deltaCol * grid.deltaCol)
 
     alpha = getAlphaX(grid)
-
     bcLeft = getBoundarySide(bc, LEFT)
     bcRight = getBoundarySide(bc, RIGHT)
+    length = getCols(grid)
 
-    concentrations::Matrix{T} = grid.concentrations[]
-    rowIndex = 1
+    concentrations::Matrix{T} = getConcentrations(grid)
 
-    numCols = max(length, 2)
+    A::Tridiagonal{T} = createCoeffMatrix(alpha, alpha_left[1, :], alpha_right[1, :], bcLeft, bcRight, length, 1, sx)
 
-    alpha_left = calcAlphaIntercell(alpha[:, 1:(numCols-2)], alpha[:, 2:(numCols-1)])
-    alpha_right = calcAlphaIntercell(alpha[:, 2:(numCols-1)], alpha[:, 3:numCols])
-
-    A::Tridiagonal{T} = createCoeffMatrix(alpha, alpha_left[rowIndex, :], alpha_right[rowIndex, :], bcLeft, bcRight, length, rowIndex, sx)
     b = concentrations[1, :]
 
-    if getType(getBoundarySide(bc, LEFT)[1]) == CONSTANT
-        b[1] += 2 * sx * alpha[1, 1] * bcLeft[1].value
+    if getType(bcLeft[1]) == CONSTANT
+        b[1] += 2 * sx * alpha[1, 1] * getValue(bcLeft[1])
     end
-    if getType(getBoundarySide(bc, RIGHT)[1]) == CONSTANT
-        b[length] += 2 * sx * alpha[1, length] * bcRight[1].value
+    if getType(bcRight[1]) == CONSTANT
+        b[end] += 2 * sx * alpha[1, end] * getValue(bcRight[1])
     end
 
-    concentrations_t1 = A \ b
-
-    concentrations[1, :] = concentrations_t1
-
-    setConcentrations!(grid, concentrations)
+    concentrations[1, :] = A \ b
 end
 
 # BTCS solution for 2D grid
@@ -220,19 +210,30 @@ end
 
 function runBTCS(grid::Grid{T}, bc::Boundary{T}, timestep::T, iterations::Int, stepCallback::Function) where {T}
     if grid.dim == 1
+        length = getCols(grid)
+        alpha = getAlphaX(grid)
+
+        alpha_left_task = Threads.@spawn calcAlphaIntercell(alpha[:, 1:(length-2)], alpha[:, 2:(length-1)])
+        alpha_right_task = Threads.@spawn calcAlphaIntercell(alpha[:, 2:(length-1)], alpha[:, 3:length])
+
+        alpha_left = fetch(alpha_left_task)
+        alpha_right = fetch(alpha_right_task)
+
         for _ in 1:(iterations)
-            BTCS_1D(grid, bc, timestep)
+            BTCS_1D(grid, bc, alpha_left, alpha_right, timestep)
 
             stepCallback()
         end
     elseif grid.dim == 2
+        rows = getRows(grid)
+        cols = getCols(grid)
         alphaX = getAlphaX(grid)
         alphaY_t = getAlphaY_t(grid)
 
-        alphaX_left_task = Threads.@spawn calcAlphaIntercell(alphaX[:, 1:(grid.cols-2)], alphaX[:, 2:(grid.cols-1)])
-        alphaX_right_task = Threads.@spawn calcAlphaIntercell(alphaX[:, 2:(grid.cols-1)], alphaX[:, 3:grid.cols])
-        alphaY_t_left_task = Threads.@spawn calcAlphaIntercell(alphaY_t[:, 1:(grid.rows-2)], alphaY_t[:, 2:(grid.rows-1)])
-        alphaY_t_right_task = Threads.@spawn calcAlphaIntercell(alphaY_t[:, 2:(grid.rows-1)], alphaY_t[:, 3:grid.rows])
+        alphaX_left_task = Threads.@spawn calcAlphaIntercell(alphaX[:, 1:(cols-2)], alphaX[:, 2:(cols-1)])
+        alphaX_right_task = Threads.@spawn calcAlphaIntercell(alphaX[:, 2:(cols-1)], alphaX[:, 3:cols])
+        alphaY_t_left_task = Threads.@spawn calcAlphaIntercell(alphaY_t[:, 1:(rows-2)], alphaY_t[:, 2:(rows-1)])
+        alphaY_t_right_task = Threads.@spawn calcAlphaIntercell(alphaY_t[:, 2:(rows-1)], alphaY_t[:, 3:rows])
 
         alphaX_left = fetch(alphaX_left_task)
         alphaX_right = fetch(alphaX_right_task)