function heat_eq(k,r)
% main code computing and ploting the solution of the heat equation as ...
% ... required in your exam problem 3. 
x=linspace(-3,3,300);
% space discretization
t=linspace(0.0001,0.1,100);
% time discretization
T=length(t);% length of the time vector to be used in the for loop
L=length(x);% length of the space vector to be used in the for loop
u=zeros(T,L);
v=zeros(T,L);
v1=zeros(T,L);
u1=zeros(T,L);
u2=zeros(T,L);
% initialization of various solutions
for i=1:T
    % for loop in time which for every moment of time t(i) will compute the
    % solution in t(i),x and store it as a row in a matrix of size T,L
      f=sqrt(4*k*t(i));
      %coefficient in the formula for the heat solution with ...
      %...diffusivity conbstant k
      f1=sqrt(4*k*r*t(i));
      %coefficient in the formula for the heat solution  with ...
      %...diffusivity conbstant k*r
    for j=1:L
        % for loop in space which for the given time t(i) computes the ...
        % ...solution in t(i), x(j) and stores it as the elemnt i,j ...
        % ... in the solution matrix of size T,L
 u(i,j)=(1/(sqrt(pi)*f)).* integral(@(s)exp(-((x(j)-s).^2)./(f^2)).*(heaviside(s+1)-heaviside(s-1)),-Inf,Inf);
 % solution computed in (t(i),x(j)) for diffusivity conbstant k
 v(i,j)=(1/(sqrt(pi)*f1)).* integral(@(s)exp(-((x(j)-s).^2)./(f1^2)).*(heaviside(s+1)-heaviside(s-1)),-Inf,Inf);
 % solution computed in (t(i),x(j)) for diffusivity conbstant k*r
    end
    u1(i,:)=exp(-t(i))*u(i,:);
    % solution computed in (t(i),x) for diffusivity conbstant k and b=1
    u2(i,:)=exp(-10*t(i))*u(i,:);
    % solution computed in (t(i),x) for diffusivity conbstant k and b=10
    v1(i,:)=exp(-t(i))*v(i,:);
    % solution computed in (t(i),x) for diffusivity conbstant k*r and b=1
end

figure(1)
% surface plots of the two solution u1 and u2 corresponding to ...
%... diffusivity conbstant k and b=1 and b=10 respectively
subplot(1,2,1)
surf(x,t,u1)
subplot(1,2,2)
surf(x,t,u2)

figure(2)
% surface plots of the two solution u1 and v1 corresponding to ...
%... b=1 and diffusivity conbstant k and k*r respectively
% use the rotation feature in the figure panel to rotate the plots for...
%...better visualization 
subplot(1,2,1)
surf(x,t,u1)
subplot(1,2,2)
surf(x,t,v1)
