f=@(x)(0);

q=@(x)(1+x.^2); a=-1;%区间左端点 b=1;%区间右端点 ya=1; yb=1; n=20; x(1)=a; y(1)=ya; x(n+1)=b; y(n+1)=yb; h=(b-a)/n;

H=zeros(n-1,n-1); g=zeros(n-1,1); for i=1:n

x(i+1)=a+h*i; end

for i=2:n-2;

H(i,i-1)=-1/h.^2; H(i,i+1)=-1/h.^2;

H(i,i)=2/h.^2+q(a+h*i); g(i,1)=f(a+h*i); end

H(1,1)=2/h.^2+q(a+h); H(1,2)=-1/h.^2; H(n-1,n-2)=-1/h.^2;

H(n-1,n-1)=2/h.^2+q(b-h); H

g(1,1)=f(a+h)+ya/h.^2; g(n-1,1)=f(b-h)+yb/h.^2; g

c=H\g; for i=2:n;

y(i)=c(i-1); end; pp=[x;y]'

plot(pp(:,1),pp(:,2),'r-');