% number of points in the exact solution plot
nplot=100;

% initialization
xplot=linspace(0,1,nplot);
uExact=sin(8*pi*xplot.*xplot);

iplot=0;
nplots=3;

% loop over the n values
for n=[10 20 50]

    % grid spacing
    dx=1/(n+1);

    % build the full matrix for the full solver
    A=(-2*diag(ones(n,1))+diag(ones(n-1,1),1)+diag(ones(n-1,1),-1))/dx/dx;

    % build the diagonals for the tridiagonal solver
    a=ones(n,1)/dx/dx;
    b=-2*ones(n,1)/dx/dx;
    c=ones(n,1)/dx/dx;

    % the grid
    x=linspace(0,1,n+2);

    %f=-64*pi*pi*sin(8*pi*x(2:n+1))
    %uExact=sin(8*pi*xplot)
    
    % the right hand side
    f=16*pi*cos(8*pi*x(2:n+1).*x(2:n+1))-256*pi*pi*x(2:n+1).*x(2:n+1).*sin(8*pi*x(2:n+1).*x(2:n+1));

    % initialize the solution vectors
    uApproxMatlab=zeros(n+2,1);
    uApproxMine=zeros(n+2,1);

    % call the matlab full matrix solver
    uApproxMatlab(2:n+1)=f/A;

    % call my tridiagonal solver
    uApproxMine(2:n+1)=tridiag(n,a,b,c,f);

    %plot the results
    iplot=iplot+1;
    subplot(nplots,1,iplot)
    plot(xplot,uExact,'-.')
    hold on
    plot(x,uApproxMatlab,'-*')
    plot(x,uApproxMine,'or')
    hold off

end

