I'm trying to test a simple 1D Poisson solver to show that a finite difference method converges with $\mathcal{O}(h^2)$ and that using a deferred correction for the input function yields a convergence with $\mathcal{O}(h^4)$.So, the equation is $ - u'' = f  $ with boundary conditions $u(0) = u(1) = 0$. The method I'm trying to use is using the discretized operator $$ A = \left[ \begin{array}{c} 2&-1&0&0&0&0 \\-1&2&-1&0&0&0 \\ 0&-1&2&-1&0&0 \\ 0&0&-1&2&-1&0 \\ 0&0&0&-1&2&-1 \\ 0&0&0&0&-1&2\end{array}\right]  $$(the example matrix is for $h = \frac{1}{5}$.) Then solve for $Au=h^2f$. I've shown that theoretically this should converge with $\mathcal{O}(h^2)$, but when I test it on Matlab, I'm getting only $\mathcal{O}(h)$ convergence.Then, I'm trying what my course instructor called deferred correction, and altering $f$ before solving. I concluded that the correction should be $f \mapsto f + \frac{h^2}{12} Af$. I've shown that this should converge with $\mathcal{O}(h^4)$, but in Matlab I still get $\mathcal{O}(h)$.Here's the Matlab script:function [u err] = threeptsolve(ureal, du2, h)% INPUT: 'ureal' is the function handle for the real solution.%        'du2' is the function handle for the second derivative of 'ureal'%        'h' is the step size% OUTPUT: 'u' is the approximated solution%         'err' is the error at each pointx = [0:h:1]';n = length(x);f = -du2(x);realu = ureal(x);A = 2 * eye(n);A = A + diag(-1*ones(n-1,1), 1) + diag(-1*ones(n-1,1), -1);A = (1/h^2) * A;% uncomment if using deferred correction% f = f + h^2/12 * A * f;u = A\f;err = (realu - u);endWhen I try this with some sample smooth functions (with appropriate boundary values), and then try again with $h/2$, I get a vector of (approximate) twos when I compute err1 ./ err2(1:2:end).Is my math wrong, or is it my code?
Testing 1D Poisson Solver
matlab;finite difference;poisson
Your formulation of $A$ assumes that $u_0$ and $u_{n+1}$ are zero, which is correct.  However, you are then including your boundaries at $u_1$ and $u_n$.  Your exact answer satisfies these later BC's but not the imposed BC's in $A$.  These discontinuities floating around cause you to lose your higher order accuracy. To get second order convergence you only need to change one line:x = [0:h:1]';tox = [h:h:1-h]';And I get a column of 4's (for 2nd order conv.) when I execute err1 ./ err2(2:2:end-1) (Note the shift by 1 index since the solutions now line up at the even indices rather than the odd).  I have not gotten 4th order yet from your deferred correction, however this solve a part of your problem.