function x = Richardson(A,b,iter,x)
  n = length(b)
  M = diag([1,1])-A
  N = diag([1,1])
  for i=1:iter
    x=M*x+N*b
end
epsilon = abs(x-inv(A)*b)
printf("Error is %f\n",max(epsilon))
endfunction

function x = Jacobi(A,b,iter,x)
  n = length(b)
  D = diag([A(1,1),A(2,2)])
  M = inv(D)*(D-A)
  N = inv(D)
  for i=1:iter
    x=M*x+N*b
end
epsilon = abs(x-inv(A)*b)
printf("Error is %f\n",max(epsilon))
endfunction

function x = GaussSeidel(A,b,iter,x)
  n = length(b)
  D = diag([A(1,1),A(2,2)])
  R = diag([0,0])
  L = diag([0,0])
  R(1,2)=A(1,2)
  L(2,1)=A(2,1)
  
  M = -inv(D+L)*R
  N = inv(D+L)
  for i=1:iter
    x=M*x+N*b
end
epsilon = abs(x-inv(A)*b)
printf("Error is %f\n",max(epsilon))
endfunction

function x = SOR(A,b,iter,x,omega)
  n = length(b)
  D = diag([A(1,1),A(2,2)])
  R = diag([0,0])
  L = diag([0,0])
  R(1,2)=A(1,2)
  L(2,1)=A(2,1)
  
  M = inv(D+omega*L)*((1-omega)*D-omega*R)
  N = omega*inv(D+omega*L)
  //printf("Maximal eigenvalue of the system matrix is %f\n",g)
  for i=1:iter
    x=M*x+N*b
end
epsilon = abs(x-inv(A)*b)
printf("Error is %f\n",max(epsilon))
endfunction

y=[21 -19]';
b=[0.3 0.3]';
A=[0.7 -0.4; -0.2 0.5];

//SOR(A,b,5,y,1.0647869)

function x=plotspec(A)
  [n,n] = size(A)
  D = diag([A(1,1),A(2,2)])
  R = diag([0,0])
  L = diag([0,0])
  R(1,2)=A(1,2)
  L(2,1)=A(2,1)
  j=1
  for (omega=0:0.001:2)  
    M = inv(D+omega*L)*((1-omega)*D-omega*R)
    v = max(real(abs(spec(M))))
    
    x(j)=v
    j=j+1
  end
  plot2d(0:0.001:2,x)

  endfunction
