clear
clc

%%%%% Inicijalizacija
n = 25;
h = 1 / (n - 1)
range = linspace(0, 4, n);
A = 2 * eye(n) - diag(ones(n - 1, 1), 1) - diag(ones(n - 1, 1), -1);
A = kron(A, eye(n)) + kron(eye(n), A);
b = zeros(n);

for i = 1 : n
  for j = 1 : n
    if i == 1 || i == n
      b(i, j) = exp(range(j)) * cos(range(i)) - exp(range(i)) * cos(range(j));
    endif
    if j == 1 || j == n
      b(i, j) = b(i, j) + exp(range(j)) * cos(range(i)) - exp(range(i)) * cos(range(j));
    endif
  endfor
endfor

b = reshape(b, n * n, 1);
eps = 1e-5;
x0 = zeros(n * n, 1);

%%%%% Plot rješenja
[X, Y] = meshgrid(range, range);
Z = exp(Y) .* cos(X) - exp(X) .* cos(Y);
surf(X, Y, Z);
%hold on

%%%%% BSO neka imamo spektralni radijus
F = diag(1 ./ diag(A)) * (diag(diag(A)) - A);
r = max(abs(eig(F)))
%%%%%

%%%%% Ubrzani jacobi
[u, res] = jacobibrzi(A, b, x0, eps, r);
size(X)
figure()
surf(X, Y, reshape(u, n, n)')