function spektralniMetoda2d_nehomogenniMaterial(k)

syms x1 x2

idx = 1;
for i=0:k
  for j=0:k-i
    Phi(idx) = x1^(i+1)*(1-x1) * x2^(j+1)*(1-x2);
    dPhidx1(idx) = diff(Phi(idx),x1);
    dPhidx2(idx) = diff(Phi(idx),x2);
    idx = idx+1;
  end
end

n = length(Phi);
A = zeros(n,n);
b = zeros(n,1);
for i=1:n
  for j=1:n
    A(i,j) = int(int(dPhidx1(j)*dPhidx1(i)+dPhidx2(j)*dPhidx2(i),x2,0,1),x1,0,1/2);
    A(i,j) = A(i,j) + 5 * int(int(dPhidx1(j)*dPhidx1(i)+dPhidx2(j)*dPhidx2(i),x2,0,1),x1,1/2,1);
  end
  b(i) = int(int(Phi(i),x2,0,1),x1,0,1);
end

alfa = A\b;
alfa'

xx1 = 0:0.1:1; xx2 = xx1;
[X1,X2] = meshgrid(xx1,xx2);
U = X1; U(:) = 0;
for j=1:n
  for k=1:length(U(:))
    Phijxk = subs(Phi(j),[x1,x2],[X1(k),X2(k)]);
    U(k) = U(k) + alfa(j)*Phijxk;
  end
end
mesh(X1,X2,U)