%Program 1 (b) 
%Heat function 2 dimension
%解線性方程是用內建函式，也可以用一些其他方法，不過因為要想其他的部分細節，
%所以這個部分就先簡單用內建
%用 (3) 的式子

%N = [100];
N = [100 150 200 250 300 350 400];
M = zeros(1,length(N));
O = zeros(1,length(N));
for n=1:7
    v = 0.1;
    ds = 1/N(n);
    dt = 1/N(n);
    r = v*dt/(ds^2);
    A = delsq(numgrid('S',N(n)-1+2)); %%內建函式存非0點，剛好是[-1 ... -1 4 -1 ... -1]
    A = (r/2)*A;
    N2 = (N(n)-1)^2;
    uo = zeros(N2,1);
    ue = zeros(N2,1);
    un = zeros(N2,1);
    %將起始值和實際值存成(N-1)^2 X 1 向量
    for i=1:N(n)-1
        for j=1:N(n)-1
            uo((i-1)*(N(n)-1)+j) = sin((i)*ds*pi)*sin((j)*ds*pi);
            ue((i-1)*(N(n)-1)+j)=exp(-2*pi*pi*v*1)*sin((i)*ds*pi)*sin((j)*ds*pi);
        end
    end
    
    I = spdiags(ones(N2,1),0,N2,N2);
    uo = (I-A)*uo;
    %疊代 N 次，一次跑dt = 1/N ，N 次 t = 1
    for k=1:N(n)
        un = (I+A)\uo;
        uo = (I-A)*un;
        %norm(un-ue,2)
    end
 
    M(n) = norm(ue-un,2); %最大誤差 2norm
    O(n) = 1/N(n);
end

%畫 最大誤差 和 1/N 的圖，若為線性應該可以表示最大誤差是 O(1/N)
%不知道我的想法有沒有錯
%下面兩個函式 plot 和 loglog 不知道有什麼差別，只是抓助教寫的作業解答
figure(1)
plot(O,M);

figure(2)
loglog(N,M);