For small deformations the equilibrium equation in the displacement formulation is [1]
(μ + λ)∇divu + μΔu
+ f = 0, f = div(I Ta) = ∇Ta,
where u = (X, Y) is displacement of a point from equilibrium position,
Ta represents the active stress in cardiac tissue [2], μ, λ are Lame parameters
(we put μ = 1 ). Using 2-nd order finite differences we obtain the approximate equation
(2 + λ)(X_{i+1,j} + X_{i-1,j} - 2X_{i,j} ) +
X_{i,j+1} + X_{i,j-1} - 2X_{i,j} + (Y_{i+1,j+1} + Y_{i-1,j-1}
- Y_{i+1,j-1} - Y_{i-1,j+1} )(1 + λ)/4 + h^{2}Fx_{i,j} = 0,
Fx_{i,j} = (Ta_{i+1,j} - Ta_{i-1,j })/2h. (*)

2D weighted Jacobi relaxation algorithm on CPU

For weighted Jacobi relaxation [3] of (X_{i,j }, Y_{i,j }) on
n×n grid we rewrite equations (*) as
2(3 + λ)X_{i,j} = (2 + λ)(X_{i+1,j} + X_{i-1,j} ) +
X_{i,j+1} + X_{i,j-1} + (Y_{i+1,j+1} + Y_{i-1,j-1}
- Y_{i+1,j-1} - Y_{i-1,j+1} )(1 + λ)/4 + h^{2}Fx_{i,j} where Ta ~ exp[-(r/a)^{2}]. 50 cross-sections X_{i, n/2} and
error log max |X_{i,j}^{new} - X_{i,j}^{old}|
(the red curve below) are ploted after every 2it iterations.
Set it = 50 to get accurate displacement values.
n
it×2
ω

You can check that relaxation is fast for small n (e.g. n = 32)
but it is very slow for fine-grids.

[1] Linear elasticity (Wiki)
[2] M.P. Nasha, A.V. Panfilov
Electromechanical model of excitable tissue to study reentrant cardiac arrhythmias Progress in Biophysics & Molecular Biology 85 (2004) 501.
[3] William L. Briggs, Van Emden Henson, Steve F. McCormick
A Multigrid Tutorial
(look at "Tutorial Slides")