Two-Point Boundary Value Problems - Numerical Methods in Engineering with MATLAB - page 315

function [c,d,e,b] = fDiffEqs(x,h,n)
% Sets up the tridiagonal coefficient matrix and the
% constant vector of the finite difference equations.
h2=h*h;
d = ones(n,1)*(-2 + 4*h2);
c = ones(n-1,1);
e = ones(n-1,1);
b = ones(n,1)*4*h2.*x;
d(1) = 1; e(1) = 0; b(1) = 0;c(n-1) = 2;

The solutionis

>> x y1

0.0000e+000 0.0000e+000

1.5708e-001 3.1417e-001

3.1416e-001 6.1284e-001

4.7124e-001 8.8203e-001

6.2832e-001 1.1107e+000

7.8540e-001 1.2917e+000

9.4248e-001 1.4228e+000

1.0996e+000 1.5064e+000

1.2566e+000 1.5500e+000

1.4137e+000 1.5645e+000

1.5708e+000 1.5642e+000

The exact solution of the problemis

y = x − sin 2 x

which yields y ( π/ 2) = π/ 2 = 1 . 57080. Thus the errorinthe numerical solutionis

about0 . 4%. More accurate results can be achievedbyincreasing n .For example, with

n = 101, we wouldget y ( π/ 2) = 1 . 57073 , which is in errorbyonly 0 . 0002%.

EXAMPLE 8.7
Solve the boundary value problem

y =− 3 yy

y (0) = 0 y (2) = 1

with the finite difference method. (This problemwas solved in Example 8.1 by the
shootingmethod.)Use n

= 11 and compare the results to the solution in Example 8.1.

Solution As the problemis nonlinear, Eqs. (8.11) must be solvedbythe Newton-
Raphsonmethod. The program listedbelow can be usedas amodel for other