function HeatDynamicImplicit()
L = 1;   % length of the space interval
n = 29;  % number of interior grid points
r = 0.2; % ratio to compute time step
%r = 20.0;% ratio to compute time step
T = 0.5;  % final time
plots = 5;% number of plots to be saved

function r = iv(x);
  r = min([ 2*x/L,2-2*x/L]')';
end%function

dx = L/(n+1);
dt = 2*r*dx^2;
x = linspace(0,L,n+2)';
initval = iv(x(2:n+1));

yplot = zeros(plots,n+2);
plotc = 1; tplot = linspace(0,T,plots);

Adiag = ones(n,1)*(1+2*r);
Aoffdiag = -ones(n-1,1)*r;

y = initval;
for t = 0:dt:T+dt;
  if min(abs(tplot-t))<dt/2
    yplot(plotc,2:n+1) = y';
    plotc = plotc +1;
  end%if
  y = trisolve(Adiag,Aoffdiag,y);
end%for
plot(x,yplot)
grid on
xlabel('position x'); ylabel('Temperature')

end%function
