function mol
u = @(t,x)0.5./sqrt(4*t+1)*(exp(3*(2*x+1)/(4*t+1))+2).*exp(-(x+2).^2/(4*t+1));
M = 40;
N = 100;
x = linspace(-5,5,M+1);
y0 = 0.5*exp(-(x-1).^2)+exp(-(x+2).^2);
[t,y] = ode45(@dmol,linspace(0,0.5,N),y0(2:end));
y = [zeros(N,1) y];

hold on
for tt=1:numel(t)
    cla
    stem(x,y(tt,:),'k')
    plot(x,u(t(tt),x),'k')
    text(0.8,0.8,['t = ' num2str(t(tt)) 's']);
    ylim([0 1])
    xlabel('x')
    ylabel('u(x,t)')
    drawnow   
end

    function dy=dmol(~,y)
        dx = 10/M;
        y0 = 0;
        for m=M:-1:1
            if m==1
                dy(m,1) = (y(m+1)-2*y(m)+y0)/(dx^2);
            elseif m==M
                dy(m,1) = 0;
            else
                dy(m,1) = (y(m+1)-2*y(m)+y(m-1))/(dx^2);
            end
        end
    end %dmol function

end %mol function