clear;
p=1000;
xl=0;xr=1;tb=0;te=1;M=p;N=p;
alf=5;bet=4;D=0.05;

%Exact solution given by (8.68) in the text
xk = [0:1/M:1];
tk = [0:1/N:1];
for i = 1:length(tk)
    for j = 1:length(xk)
        u(i,j) = (2*D*bet*pi*exp(-D*pi^2*tk(i))*sin(pi*xk(j)))...
            /(alf+bet*exp(-D*pi^2*tk(i))*cos(pi*xk(j)));
    end
end

surf(xk,tk,u'); % 3-D plot of solution w
hold on;
%shading interp %Uncomment these to make the surface look pretty
%material dull
title('Exact solution: \alpha = 5, \beta = 4, D = 0.05','FontSize',15)
xlabel('x position (units)','FontSize',15)
ylabel('time (units)','FontSize',15)
zlabel('u','FontSize',15)
hold off;
shading interp
material dull
