%clear;
close all;
h=10;k=30;md=7300;cp=800;t0=1000;tw=0;tsp=600;le=0.1;%tw為環境溫度
n_space = input('n_space?   ');
n_space=n_space+1;
n_time = input('n_time?   ');
alfa=tsp/n_time/(le*le/n_space/n_space)*k/md/cp;
alfa1=1+2*alfa;

beta=(le/n_space)*h/k;

x(1,1)=tw*beta;
x(n_space,1)=tw*beta;
x(2:(n_space-1),1)=t0;
b=x;
%矩陣計算:
D(n_space,n_space)=1+beta;
D(1,1)=1+beta;
L(n_space,n_space)=0;
L(n_space,n_space-1)=1;
U(n_space,n_space)=0;
U(1,2)=1;

for ii=2:(n_space-1)
    D(ii,ii)=alfa1;
    L(ii,ii-1)=alfa;
    U(ii,ii+1)=alfa;   
end
J=pinv(D)*(L+U);pj=vrho(J);

w=2/(1+sqrt(1-(pj*pj)));
clear J;
MM=pinv(D-w*L);
loopmatrix=MM*((1-w)*D+w*U);     %%%%%%%%%%%迭代矩陣
MM=w*MM;
    
    for tt=1:n_time%對時分的循環
          zzz=1;stop=1;
         while stop>0.01%某時間片一次SOR循環 循環控制條件 有待改進~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
              
              x1=loopmatrix*x+MM*b;
              stop=norm(x-x1,2);
              x=x1; 
              disp(['時間步= ' , int2str(tt), ';第幾次SOR： ' ,int2str(zzz)]); stop
              zzz=zzz+1;
         end
         count(tt)=zzz;
          if tt<n_time
              b=x;
              b(1,1)=0;
              b(n_space,1)=0;
              
          end
    
    end
    sor_average_times=mean(count)
    w
    pj
    for i=1:n_space length_dis(i)=(le/(n_space-1))*(i-1);end
    x2=x;length_dis2=length_dis;
    figure;
plot(length_dis,x,'k');
title (['\fontsize{12}\bf \rm','溫度場曲線圖']);
xlabel ('棒子空間，以長度計量'); ylabel ('溫度');