function [T,timeV,spaceV]=nma_PDE_parabolic_explicit_rod_with_rate_BC(length,leftTRate,rightTRate,rodT,deltaT,nSteps,deltaX,k)
%solve parabolic PDE using explicit method for changing boundary conditions
%function [T,timeV,spaceV]=nma_PDE_parabolic_explicit_rod_with_rate_BC(length,leftTRate,rightTRate,rodT,deltaT,deltaX,k)
%
%Solves the parabolic PDE using explicit method for given Boundary
%expressed as rate (dT/dt). i.e. B.C. is not fixed values.
%
%
%INPUT:
% length: length of rod in some length units
% leftTRate: rate of change of independent variable (V or T for example) at
% left end of rod w.r.t time
% rightTRate: rate of change of independent variable (V or T for example) at
% right end of rod w.r.t. time
% 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);
%set initial conditions;
T(1,:)=rodT;
timeV=zeros(nSteps,1);
spaceV=zeros(nPoints,1);
for(i=1:nPoints)
spaceV(i)=(i-1)*deltaX;
end
lambda=k*deltaT/(deltaX)^2;
for(t=1:nSteps-1)
for(i=2:nPoints-1)
%do left edge
T(t+1,1)=T(t,1)+leftTRate*deltaT;
T(t+1,end)=T(t,end)+rightTRate*deltaT;
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