%Program 2 (b) Code
%作法和2 (a) 類似，差別在於 |U|^2 值，這部分和Program 1 (c)的問題一樣，
%不確定|U|^2是要怎麼存取，有兩種想法
%第一種：切割完分割點後，對每個分割點的 U(x_i,y_j) 各別平方，然後存成 (N-1)^2 X 1
%        的向量
%第二種：取 Norm^2
%這裡我是用第一種方法，且取 U(x_i,y_j)=1，這樣可能才會符合那個疊積分限制式
%如果我的想法不對的話，應該只需要重新設 |U|^2 值

N = 200;
ds = 1/N;
r = 1/ds^2;
R = delsq(numgrid('S',N-1+2)); %內建函式存非0點，剛好是[-1 ... -1 4 -1 ... -1]
R = r*R;
N2 = (N-1)^2;
uo = zeros(N2,1);
uo2 = zeros(N2,1);
un = zeros(N2,1);
for i=1:N-1
    for j=1:N-1
        %uo((i-1)*(N-1)+j) = sin((i)*ds*pi)*sin((j)*ds*pi);
        uo((i-1)*(N-1)+j) = 1;
        uo2((i-1)*(N-1)+j) = uo((i-1)*(N-1)+j)^2;
    end
end
v = zeros(N2,1);
alpha = 1;
for i=1:N-1
    for j=1:N-1
        v((i-1)*(N-1)+j) = alpha*sin((i)*ds*pi)*sin((j)*ds*pi);
    end
end

V = spdiags(v+uo2,0,N2,N2);
A = R+V;
w = zeros(N2,1);
eo = 0;
k = 1;
tol = 10^-5;
while k<=1000
    w = A\uo;
    un = w/norm(w,2);
    en = un'*A*un;
    if abs(eo-en)<tol
        en
        break
    end
    k = k+1;
    eo = en;
    uo = un;
end