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

function JacobiP(A)
  D = diag([A(1,1),A(2,2)])
  B = inv(D)*(D-A)
  PL = inv(D)
  c1=cond(D-A)
  c2=cond(B)
  disp(D)
  disp(PL)
  disp(A)
  disp(PL*A)
  disp(cond(PL*A))
  disp(cond(A))
endfunction

function GaussSeidelP(A,b,iter,x)
  D = diag([A(1,1),A(2,2)])
  L = diag([0,0])
  L(2,1)=A(2,1)
  PL = inv(D+L)
  disp(D)
  disp(PL)
  disp(A)
  disp(PL*A)
  disp(cond(PL*A))
  disp(cond(A))
endfunction

function SORP(A,omega)
  D = diag([A(1,1),A(2,2)])
  L = diag([0,0])
  L(2,1)=A(2,1)
  PL = omega*inv(D+omega*L)
  disp(D)
  disp(PL)
  disp(A)
  disp(PL*A)
  disp(cond(PL*A))
  disp(cond(A))  
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
