J=16
h = 1/J
for i=1:J+1 
int(i) =(i-1)*h ;
end
for j=1:J-1
uexact(j) = (8/pi^3)*(exp(-pi^2)*sin(pi*j*h) - exp(-9*pi^2)*sin(3*pi*j*h)/27);
end
% Define the exact solution including the boundary values
uexactp = [0,uexact,0];
plot(int,uexactp)
% Note that if the approximate solution un is defined at the mesh points
% j*h, j=1,...J-1, then defining
% unp=[0,un,0]
% both unp and uexactp can be plotted on the same graph by the command
% plot(int,unp,int,uexactp)