%Program 2 (a) Code
%解線性方程是用內建函式，也可以用一些其他方法，不過因為要想其他的部分細節，
%所以這個部分就先簡單用內建
N = 100;
ds = 1/N;
r = 1/ds^2;
R = delsq(numgrid('S',N-1+2)); %%內建函式存非0點，剛好是[-1 ... -1 4 -1 ... -1]
R = r*R; %係數矩陣Uxx+Uyy
N2 = (N-1)^2;
uo = zeros(N2,1);
un = zeros(N2,1);

%V(x,y)存成一行向量
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
%起始值U^(0)，不知道設多少，不過用以下兩種起始值結果好像一樣。
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;
    end
end
%將V(x,y)變成方陣
V = spdiags(v,0,N2,N2);
A = R+V;  %左式U^(1) 的係數矩陣
%作(11)公式的迴圈
w = zeros(N2,1);
lamda0 = 0; % 假定起始Lamda是0，不是很確定要設多少
k = 1;
tol = 10^-5;
while k<=1000
    w = A\uo;
    un = w/norm(w,2);
    lamda1 = un'*A*un;
    %前後Lamda值小於誤差就印出lamda1，這和lamda0有關，不過跑出來好像沒影響
    if abs(lamda0-lamda1)<tol
        lamda1
        break
    end
    k = k+1;
    lamda0 = lamda1;
    uo = un;
end