spectral relaxation method for ode

I'm trying to create a function that performs the SOR method. how to give an input matrix which is given below..

A1=(1+1/beta)D^2+diag(f_r+g_r)D-(2h_r+M)I

B1=2P_r^2

initial guess

f_0=1-exp(-eta)

p_0=exp(-eta)

g_0=alpha(1-exp(-eta))

h_0=alpha(exp(-eta))

where beta=0.8,M=0.5