function analytical_solution_of_u(N,nn)
% N represents the first N terms in the summation of the series.
% nn is the number of the mesh points 
%g=e^(-200*(s - 1/2).^2)
%C_S=cn*sin(n*π*x/L)   E=(-(n^2π^2/L^2)^(α/2)*t)
%uexact=∑^{x=N}_{x=1} C_S*sin(n*π*x)*E
Alp=[1.1,1.3,1.5,1.8]; % α 
T=0.5; %end time
Eps = 0.01; %u_t = -Eps(-Δ)^(Alp/2)u
C_S = @(n)2*integral(@(s)exp(-200*(s - 1/2).^2).*sin(n*pi.*s),0,1);
E = @(alp,t,eps,n)exp(-eps*(n^2*pi^2).^(alp/2)*t);
for i=1:N
    cn(i) = C_S(i);
    Exp(1:4,i) = E(Alp,T,Eps,i);
end
x = linspace(0,1,nn);
for j = 1:length(Alp)
    for i=1:size(x,2)
        s(i,:) = sin((1:N)*pi*x(i));
        ue(j,i) = sum(s(i,:).*cn.*Exp(j,:));
    end
    plot(x,ue(j,:),'LineWidth',2.15);
    hold on
end
legend('\alpha=1.1','\alpha=1.3','\alpha=1.5','\alpha=1.8');
set(gca,'FontSize',23)
xlabel('$x$ ','interpreter','latex','FontSize',33)
ylabel('$u$ ','interpreter','latex','FontSize',33,'Rotation',0)
axis([0 1 0 1]);
set(gcf,'position',[300,100,900,700])
