%
% Temperature distribution in a 'furnace supporting column'
%  n = number of gridpoints in x-direction
%  ome = relaxation parameter,  1<=ome<2  
%
% Laplace(u)=0  0<x<1/2, 0<y<1
% u(x,y)=500  if x=0 or y=0
% u_x(x,y)=0 if x=1/2
% u_y(x,y)+10u(x,y)=3000 y=1
%
function furnace(n,ome)
m=2*n+1;
h=1/(m+1);
h2=h*h;
xx=0:h:0.5;
yy=0:h:1;
[X,Y]=meshgrid(xx,yy);
%f=zeros(m+2,n+2);
u=500*ones(m+2,n+2);
itmax=1000;
ome1=1-ome;
for iter=1:itmax
    uold = u;
    for i=2:m+1
        %
        % Inner nodes
        %
        for j=2:n+1
            uu = 0.25*(u(i+1,j)+u(i-1,j)+u(i,j+1)+u(i,j-1));
            u(i,j) = ome1*u(i,j) + ome*uu;
        end
        %
        % Neumann nodes at x=1/2
        %
        uu = 0.25*(2*u(i,n+1)+u(i-1,n+2)+u(i+1,n+2));
        u(i,n+2) = ome1*u(i,n+2) + ome*uu;
    end
    %
    % Robin nodes at y=1:
    %
    for j=2:n+1
        uu=(2*h*3000+u(m+2,j+1)+u(m+2,j-1)+2*u(m+1,j))/(4+10*2*h);
        u(m+2,j) = ome1*u(m+2,j) + ome*uu;
    end
    %
    % Upper right corner
    %
    uu=(2*h*3000+2*u(m+2,n+1)+2*u(m+1,n+2))/(4+10*2*h);
    u(m+2,n+2) = ome1*u(m+2,n+2) + ome*uu;
    %
    % Check convergence:
    %
    if norm(u-uold,'inf') < 1e-8*norm(u,'inf')
        break;
    end
end
fprintf(1,'S.O.R. iterations needed = %d ...\n',iter);

subplot(1,2,1)
[c,h]=contourf(X,Y,u,10); axis equal; title('temperature distribution');
set(h,'EdgeColor','none'); 
%
p=get(gca,'position'); % save position
colorbar;
set(gca,'position',p); % restore position
subplot(1,2,2)
plot(0.5*(0:n+1)/(n+1),u(m+2,:)); title('surface temperature at y=1');