%Program 1 (c) Code
%這裡我不是很確定 U^3 要怎麼存取，我有兩種想法
%第一種：切割完分割點後，對每個分割點的U(x_i,y_j)各別平方，然後存成(N-1)^2 X 1
%        的向量
%第二種：先取計算 (N-1) X (N-1) 的 U^3 方陣，再排成 (N-1)^2 X 1 的向量
%這裡我是用第一種方法
%如果我的想法不對的話，應該只需要重新設 U^3 值

N = 100;
ds = 1/N;
r = 1/(ds^2);
A=delsq(numgrid('S',N-1+2)); %內建函式存非0點，剛好是[-1 ... -1 4 -1 ... -1]
A=-r*A;
N2 = (N-1)^2;
uo = zeros(N2,1);
uo2 = zeros(N2,1);
uo3 = 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)/(2*pi*pi);
        uo2((i-1)*(N-1)+j) = uo((i-1)*(N-1)+j)^2;
        uo3((i-1)*(N-1)+j) = uo((i-1)*(N-1)+j)^3+sin((i)*ds*pi)*sin((j)*ds*pi);
    end
end
I = spdiags(3*uo2,0,N2,N2);
%運用 Newton's Method ，根據課本Ch10的定理，應該能保證收斂是二階
k = 1;
while k<=1000
    F = A*uo-uo3;
    J = A+I;
    un = -J\F;
    uo = un+uo;
    if norm(un,2)<10^-5
        uo(N2,1)  %不清楚要印什麼出來比較
        break
    end
    k = k+1;
    for i=1:N-1
        for j=1:N-1
            uo2((i-1)*(N-1)+j) = uo((i-1)*(N-1)+j)^2;
            uo3((i-1)*(N-1)+j) = uo((i-1)*(N-1)+j)^3+sin((i)*ds*pi)*sin((j)*ds*pi);
        end
    end
end