function [T,timeV,spaceV]=nma_PDE_parabolic_explicit_rod(length,leftT,rightT,rodT,deltaT,nSteps,deltaX,k)
%solves parabolic PDE using explicit method

%function [T,timeV,spaceV]=nma_PDE_parabolic_explicit_rod(length,leftT,rightT,rodT,deltaT,deltaX,k)
%
%Solves the parabolic PDE using explicit method for given Boundary
%fixed conditions (do not change with time).
%
%INPUT:
%  length: length of rod in some length units
%  leftT:  Value of independent variable (V or T for example) at 
%          left end of rod
%  rightT: Value of independent variable (V or T for example) at 
%         right end of rod
%  rodT:  Value of independent variable (V or T for example) at 
%         rod interior points at t=0 (initial)
% deltaT:  time step size
% nSteps:  number of time steps to run for.
% deltaX:  step size in distance along rod
% k: material coefficent. depends on problem (heat, Voltage, etc...)
%
% OUTPUT:
%     T : a matrix that represent the solution. Each row is the solution
%         at different time. i.e. first row is the solution at t=0;
%         second row is the solution after one time step, etc...
%         first column is solution at left edge, second column is solution
%         at one distance step to the right of the left edge, etc... 
%         until last column is reached, which is solution at right edge.
%    timeV : a vector that has the time values corresponding to when each
%        row in T was generated.
%    spaceV : a vector that has the space values corresponding to the points 
%             on the rod. These are the points where each column in T is
%             located.
%
%Author: Nasser Abbasi

nIntervals=length/deltaX;
nPoints=nIntervals+1;
T=zeros(nSteps,nPoints);
timeV=zeros(nSteps,1);
spaceV=zeros(nPoints,1);
for(i=1:nPoints)
 spaceV(i)=(i-1)*deltaX;    
end

%set initial conditions;
T(1,2:end-1)=rodT;

%set boundary conditions
T(:,1)=leftT;
T(:,end)=rightT;

lambda=k*deltaT/(deltaX)^2;

for(t=1:nSteps-1)
    for(i=2:nPoints-1)
        T(t+1,i)=T(t,i)+lambda*(T(t,i+1)-2*T(t,i)+T(t,i-1));        
        timeV(t+1)=timeV(t)+deltaT;         
    end
end