function outvar = Steady_State_Heat(a,b,BC_a,BC_b,n)
% STEADY_STATE_HEAT: equation in 1 dimension
% from Frederic Gibou's notes

x = linspace(a,b,n);
dx = (b-a)/(n-1);
A = zeros(n,n);
A(1,1) = 1; 
A(n,n) = 1;
rhs(1) = BC_a; 
rhs(n) = BC_b;

for i = 2:n-1
    A(i,i)   =  2/dx/dx; 
    A(i,i-1) = -1/dx/dx; 
    A(i,i+1) = -1/dx/dx;
    rhs(i)   = 0;    % no internal source of heat
    & rhs(i) = -(x(i)+3) * exp(x(i)); % a nonzero "heat source" term
end

outvar = A\rhs';
exact = nan(size(x));
% exact = (x+1).*exp(x); % the exact solution with the nonzero source term

plot(x,exact,'r',x,outvar,'bo');
legend('Exact','Numerical');
xlabel('x');ylabel('Temperature');
title('Numerical Solution of the Temperature Distribution');
end
