k = 1;
x = linspace(0,10,50);
u0 = exp(-(x-5).^2);
ub = [0 0];

dx = x(2)-x(1);

for steps=223:230
    t = linspace(0,5,steps);
    dt = t(2)-t(1);
    C = k*dt/dx/dx;
    u = zeros(length(t),length(x));
    u(1,:) = u0;
    u(2:end,1) = ub(1);
    u(2:end,end) = ub(2);
    for n=1:length(t)-1
        u(n+1,2:end-1) = (1-2*C)*u(n,2:end-1) + C*(u(n,3:end)+u(n,1:end-2));
    end
    surf(x,t,u)
    zlim([-0.5 1.5])
    title(['C = k \Deltat/\Deltax = ' num2str(C)],'FontSize',14,'FontWeight','bold')
    xlabel('x','FontSize',14,'FontWeight','bold')
    ylabel('t','FontSize',14,'FontWeight','bold')
    set(gca,'FontSize',14,'FontWeight','bold')
    shading interp
    drawnow
    pause(1)
end