% Labuppgift 1, värmeledning i platta

m=20;
h=1/(m+1);
n=2*m+1;
x=(0:h:2);
y=(0:h:1);
[X,Y]=meshgrid(x,y);        % Skapar koordinatsytemet på plattan

% Skapa matris
e=ones(m*n,1);
e1=ones(m*n,1);
for d=1:41
e1(m*d+1)=0;                % Förskjut radvärdena
end
e1=e1(1:m*n);               % Trimma längden

e2=ones(m*n,1);
for d=1:41
e2(m*d)=0;
end

A=spdiags([-e -e2 4*e -e1 -e],[-m -1 0 1 m],m*n,m*n);

% Skapa kolonvektor med randtemperaturer
R=zeros([1,820]);
R(1)=120;
R(2:19)=100;
for i=1:40
    R(20*i)=60;
    R(20*i+1)=20;
end
R(20)=160;
R(800)=60;
R(801)=100;
R(802:819)=80;
R(820)=140;

T=A\R';

Z=reshape(T,m,n);

% Ordna randtemperaturer
a=ones(20,1); a=80.*a;
b=ones(42,1); b=20.*b;
c=ones(21,1); c=100.*c;
d=ones(43,1); d=60.*d;

Z=[Z a];                       % Tillägg av randtemp på matrisen
Z=[b'; Z];
Z=[c Z];
Z=[Z; d'];

v=[20 30 40 50 60 70 80 90];   % Vektor med isotermer

figure(1)
surf(X,Y,Z)

figure(2)
[Z,s]=contour(Z,v);
clabel(Z,s,v)

figure(3)
spy(A,'o')