findiff.m

国外计算机科学教材系列 numerial method using matlab (third edition)

function F=findiff(p,q,r,a,b,alpha,beta,N)
%Input - p, q, and r are the coefficient functions of (1)
%        input as strings; 'p',  'q', 'r'
%      - a and b are the left and right endpoints
%      - alpha =x(a) and beta=x(b)
%      - N is the number of steps
%Output - F=[T' X']:wherer T' is the 1xN vector of abscissas and
%         X' is the 1xN vector of ordinates.

%Initialize vectors and h

T=zeros(1,N+1);
X=zeros(1,N-1);
Va=zeros(1,N-2);
Vb=zeros(1,N-1);
Vc=zeros(1,N-2);
Vd=zeros(1,N-1);
h=(b-a)/N;

%Calculate the constant vector B in AX=B

Vt=a+h:h:a+h*(N-1);
Vb=-h^2*feval(r,Vt);
Vb(1)=Vb(1)+(1+h/2*feval(p,Vt(1)))*alpha;
Vb(N-1)=Vb(N-1)+(1-h/2*feval(p,Vt(N-1)))*beta;

%Calculate the main diagonal of A in AX=B

Vd=2+h^2*feval(q,Vt);

%Calculate the super diagonal of A in AX=B

Vta=Vt(1,2:N-1);
Va=-1-h/2*feval(p,Vta);

%Calculate the sub diagonal of A in AX=B

Vtc=Vt(1,1:N-2);
Vc=-1+h/2*feval(p,Vtc);

%Solve AX=B using trisys

X=trisys(Va,Vd,Vc,Vb);
T=[a,Vt,b];
X=[alpha,X,beta];
F=[T' X'];