function [S,u] = heat_eq(k,l)
x=linspace(-1,1,1000);
x1=linspace(-5,5,100);
t=linspace(0,0.01,100);
t1=linspace(0,5,100);
T=length(t);
L1=length(x1);
L=length(x);
%Q=zeros(T,L);
S=zeros(T,L);
u=zeros(T,L1);
I=zeros(1,T);
for i=1:T
    f=sqrt(4*k*t(i));
      f1=sqrt(4*k*t1(i));
      %f2=sqrt(4*k*(t1(i)+0.01))
     S(i,:)=1/(sqrt(pi)*f).*exp(-(x.^2)./(f^2));
     %I(i)=integral(@(s)1/(sqrt(pi)*f2).*exp(-(s.^2)./(f2^2)),-Inf,Inf);
    for j=1:L1
%Q(i,j)=1/2+(1/sqrt(pi) ) * quad(@(s)exp(-s.^2),0,x(j)/f);
 u(i,j)=1/(sqrt(pi)*f1).* integral(@(s)exp(-((x1(j)-s).^2)./(f1^2)).*(exp(-s).*heaviside(s)),-Inf,Inf);
    end
    
end
%exp(-s).*heaviside(s)
%heaviside(s+l/10)-heaviside(s-l/10)+heaviside(s-2+l/10)-heaviside(s-2-l/10)
figure(1)
surf(x,t,S)

figure(2)
surf(x1,t1,u)
