
tvals = linspace(0,4,61);
xvals = linspace(-5,5,21);

init = 20 + 5*sign(xvals);  

%Initial Condition
% init = 0*xvals;
% init(11:15)=xvals(11:15);
% init(16:21)=5-xvals(16:21);

%Solving the Partial Difference Equation
uvals = heatdiscrete(1, xvals, tvals, init, [15 25]);
surf(xvals, tvals, uvals)
xlabel x; ylabel t; zlabel u  
pause

% FRame by Frame movie
Nframes=61;
heatmovie=moviein(Nframes);

snap0 = uvals(1,:);
for j=1:Nframes; 
    snap1 = uvals(j,:);
    % disp('Computed at time t='),(j-1)*dt,
    plot(xvals,snap0,xvals,snap1,'o-'),
    title(strcat('s=',num2str(s))), hold off
    movie(:,j)=getframe;
end