function  heat(nt,nx)
% a typical run command is    heat(200,10)

% copy this file to Matlab editor and check for correct file transfer !
% make sure to copy the entire file, not just what you see!

% this program uses the most basic elements available in all programming languages
% as such it does not take advantage of MatLab parallel capabilities 

% Basic solver for the heat equation     U_t= alpha U_xx
% with Dirichlet boundary conditions on  0<t<T  0<x<L
% U=U(t,x)  temperature
% U(0,x)=initial(x)
% U(t,0)=left(t)
% U(t,L)=right(t)
% nt number of intervals for time,  nt+1 time values
% nx number of intervals for space, nx+1 space values
% ht = T/nt  time-step
% hx = L/nx  space-step
% it time index
% ix space index
% s = alpha*ht/hx^2,    s  should be smaller than 1/2 for stability
% u(it,ix) is the temperature matrix corresponding to the time-space point    t=(it-1)*ht, x=(ix-1)*hx 
% u(it,ix) ~ U((it-1)*ht, (ix-1)*hx)

% Three basic but troubling issues to remember in mathematical programming
% 1- Matlab starts its counters at 1
% 2- If you divide a line to  10 segments you have 11 end points
% 3- When you discretize a continuous equation you have to convert between cell numbers and physical measurements

%input parameters
global L T alpha  % these variables will be available to subroutines below or command window
                  % care to be taken not to reset them on command window or other functions

alpha=1;
L=1;
T=1;

%derived parameters
ht=T/nt;
hx=L/nx;
s=alpha*ht/hx^2;

% reserving memory while initializing the matrix to zero
u=zeros(nt+1,nx+1); 

% checking that input parameters lead to a stable solution
if(s>=1/2)
    nx
    hx
    nt
    ht
    s
    Nt=ceil(2*alpha*T/hx^2); %  
    fprintf('unstable, increase nt to at least %i , stopping \n',Nt)
    return
end

% setting the initial condition
for ix=1:nx+1
    x=(ix-1)*hx;
    u(1,ix)=initial(x);
end

% setting the Dirichlet boundary conditions
for it=1:nt+1
    t=(it-1)*ht;
    u(it,1)=left(t);
    u(it,nx+1)=right(t);
end

% main updating step
for it=2:nt+1
    for ix=2:nx
        u(it,ix)=u(it-1,ix)+s*(u(it-1,ix-1)-2*u(it-1,ix)+u(it-1,ix+1));
    end
end

% preparing for graphing
% you can click on rotation icon on graph window to use mouse to rotate the graph
xp=linspace(0,L,nx+1);
tp=linspace(0,T,nt+1);
surf(xp,tp,u)

end % end of function heat

% definining the initial condition function
function y=initial(x)
global L T alpha
y=x*(L-x);
end

% definining the left boundary condition function
function y=left(t)
global L T alpha
y=0;
% c= alpha*T/L^2
% y= sin(3*c^2*t/T)  % an alternate condition
end

% definining the left boundary condition function
function y=right(t)
global L T alpha
y=0;
% c= alpha*T/L^2
% y= sin(2*c^2*t/T)  % an alternate condition
end