%Temp distribution in a cylinder
h=1;
r=[0:h:10];
n=length(r);
%T(1)=120;
T(n)=60;
R=r(1:(n-1));
nn=length(R);

A=zeros(nn,nn);
B=-(ones(nn,1));

for k=1:nn
A(k,k)=-2/h^2;

if k==1
%do nothing
else
A(k,k-1)=(1/h^2-1/(2*h*R(k)));
end

if k==nn
B(k,1)=-1-(1/h^2+1/(2*h*R(k)))*T(n);
%A(k,k-1)=(1/h^2-1/(2*h*R(k)))+(1/h^2+1/(2*h*R(k)));
elseif k==1
A(k,k+1)=2/h^2;
else
A(k,k+1)=(1/h^2+1/(2*h*R(k)));
end

end