function FiniteDifferenceMethod(nnodes)
% Simple Matlab program to demonstrated the finite difference
% method with the heat transfer example problem used in lecture.
% The differential equation for this problem is
% d^2T/dx^2 + h(Ta - T) = 0
% with x = 0 to 10, Ta = 20, and boundary conditions
% T(0) = 40 and T(10) = 200.
% For this approach, there is no need to convert the second order
% ODE into two first order ODEs. Instead, we will discretize the
% ODE by replacing the derivative terms with central difference
% equations:
% d^2T/dx^2 = (T(i+1) - 2*T(i) + T(i-1))/deltax^2
% Substituting this approximation into the original ODE produces
% (T(i+1) - 2*T(i) + T(i-1))/deltax^2 + h*(Ta - T(i)) = 0
% Collecting like terms of T yields
% -T(i-1) + (2+h*deltax^2)*T(i) - T(i+1) = h*deltax^2*Ta
% We now discretize the rod into a number of nodes in the x direction
% from x = 0 to x = 10 and apply the discretization equation above
% to each node.
% Note that the first and last node are not used explicitly, since we
% don't have T(i-1) and T(i+1) for those nodes). Also, we need to treat
% the second and next to last node differently. For the second node,
% we know that T(i-1) = 40, so T(i-1) is not an unknown and gets moved
% to the right side of the equation. Similarly for the next to last node,
% we know that T(i+1) = 200, so T(i+1) is not an unknown and gets moved
% to the right side of the equation.

% Specify known problem constants
h = 0.01;
Ta = 20;
T0 = 40;
T10 = 200;

% Discretize the problem based on the specified number of nodes
n = nnodes-2;
deltax = 10/(nnodes-1)
A = zeros(n,n);
b = zeros(n,1);

% Form linear system of equations to be solved
for i = 1:n
    if i == 1
        A(1,1) = 2+h*deltax^2;
        A(1,2) = -1;
        b(1,1) = h*deltax^2*Ta+T0;
    elseif i == n
        A(n,n-1) = -1;
        A(n,n) = 2+h*deltax^2;
        b(n,1) = h*deltax^2*Ta+T10;
    else
        A(i,i-1) = -1;
        A(i,i) = 2+h*deltax^2;
        A(i,i+1) = -1;
        b(i,1) = h*deltax^2*Ta;
    end
end

A
b

% Solve linear system for unknown nodal temperatures
T = zeros(nnodes,1);
T(1,1) = T0;
T(nnodes,1) = T10;
T(2:nnodes-1,1) = A\b

x = linspace(0,10,nnodes)';
plot(x,T)
xlabel('x')
ylabel('T')