I'm learning finite difference method for solving boundary value problems. So, I want to solve the following bvp using finite difference method, especially I'm using forward difference. My intention is to show the difference accuracy of forward difference and centered difference. The bvp is $$ \begin{array}{*{20}{c}} {y''-5y'+6y=2{{e}^{t}},} & {y'\left( 0 \right)=1,y\left( 2 \right)={{e}^{2}}} \end{array}$$
Using forward difference I get: $$ \begin{array}{l}y''-5y'+6y=2{{e}^{t}}\\\frac{1}{{{{h}^{2}}}}\left( {{{y}_{{n+2}}}-2{{y}_{{n+1}}}+{{y}_{n}}} \right)-\frac{5}{h}\left( {{{y}_{{n+1}}}-{{y}_{n}}} \right)+6{{y}_{n}}=2{{e}^{t}}\\\left( {\frac{1}{{{{h}^{2}}}}+\frac{5}{h}+6} \right){{y}_{n}}+\left( {\frac{{-2}}{{{{h}^{2}}}}-\frac{5}{h}} \right){{y}_{{n+1}}}+\left( {\frac{1}{{{{h}^{2}}}}} \right){{y}_{{n+2}}}=2{{e}^{t}}\end{array}$$
So, I let ${{a}_{1}}=\frac{{1+5h+6{{h}^{2}}}}{{{{h}^{2}}}};{{a}_{2}}=\frac{{-2-5h}}{{{{h}^{2}}}};{{a}_{3}}=\frac{1}{{{{h}^{2}}}}$
so, for $n=2$ until $n=N-2$, the equation will be $$ {{a}_{1}}{{y}_{n}}+{{a}_{2}}{{y}_{{n+1}}}+{{a}_{3}}{{y}_{{n+2}}}=2{{e}^{t}}$$
For $n=N-1$, I apply backward difference method, since $y_{N+1}$ is unknown if I use forward difference, applying backward difference I get: $$ \begin{array}{l}\frac{1}{{{{h}^{2}}}}\left( {{{y}_{{n-1}}}-2{{y}_{{n-2}}}+{{y}_{{n-3}}}} \right)-\frac{5}{h}\left( {{{y}_{{n-1}}}-{{y}_{{n-2}}}} \right)+6{{y}_{{n-1}}}=2{{e}^{t}}\\\left( {\frac{1}{{{{h}^{2}}}}} \right){{y}_{{n-3}}}+\left( {\frac{{-2}}{{{{h}^{2}}}}+\frac{5}{h}} \right){{y}_{{n-2}}}+\left( {\frac{1}{{{{h}^{2}}}}-\frac{5}{h}+6} \right){{y}_{{n-1}}}=2{{e}^{t}}\end{array}$$
and in boundary : for $n=1$ $$ \begin{array}{l}{{a}_{1}}{{y}_{n}}+{{a}_{2}}\left( {{{y}_{n}}+h} \right)+{{a}_{3}}{{y}_{{n+2}}}=2{{e}^{t}}\\\left( {{{a}_{1}}+{{a}_{2}}} \right){{y}_{n}}+{{a}_{3}}{{y}_{{n+2}}}=2{{e}^{t}}-{{a}_{2}}h\end{array}$$ for $n=N+1$, the matrix element will be 1.
So, applying in matlab, I have:
function y=forward(h)
%y"-5y'+6y=2*exp(t);y'(0)=1;y(2)=exp(2)
x=0:h:2; N=length(x); t=linspace(0,2,N);
%exact solution ye=exp(x);
A=zeros(N,N); b=zeros(N,1);
%value a1,a2,a3 for n=2 until n=N-2
a1=(1/(h^2))+(5/h)+6;
a2=(-2/(h^2))-(5/h); a3=1/(h^2);
%boundary condition A(1,1)=a1+a2; A(1,3)=a3; A(N,N)=1;
%for n=N-1, using backward difference A(N-1,N-3)=1/h^2;
A(N-1,N-2)=(-2/h^2)+(5/h);
A(N-1,N-1)=(1/h^2)-(5/h)+6;
b(1)=2*exp(x(1))-a2*h;
b(N-1)=2*exp(x(N-1));
b(N)=exp(2);
for n=1:N-2 A(n,n)=a1;
A(n,n+1)=a2; A(n,n+2)=a3; b(n)=2*exp(t(n));end
y=A\b;
end
After that, I make a script to make a graph for various value of h, and here is the result I get:
Even, after using h=0.001, the result is not really good. I have no idea where I did wrong, am I missing something? Thanks a lot :)