bvp1.m

来自「用matlab写的一些数值算法」· M 代码 · 共 22 行

function  [x,y] = bvp1(q,f,bv,ab,N)
% Solve the boundary value problem  y'' - q(x)y = f(x) ,
% y(a) = bv(1), y(b) = bv(2);   a=ab(1), b=ab(2).
% Grid with step length  h = (b-a)/N.
% q(x)  is assumed to be nonnegative.

% Version 11.12.2003.  INCBOX

a = ab(1);  b = ab(2);  
h = (b - a)/N;  
x = linspace(a,b,N+1)';   % grid.  Column vector
% Function values at interior points
qv = feval(q,x(2:N));  fv = feval(f,x(2:N));
% Tridiagonal matrix in sparse format
A = spdiags([ones(N-1,1) -(2+h^2*qv) ones(N-1,1)], ...
  -1:1, N-1,N-1);
% Right hand side
b = h^2*fv;
b(1) = b(1)-bv(1);  b(end) = b(end)-bv(2);
% Solve the system and append the boundary values
y = [bv(1); A\b; bv(2)];