function [S] = swanju( x,f,f1a,f1b )
% swanju 第一种边界条件的三弯矩方法:
% 此处显示详细说明
% x 为自变量x
% f 为x对应的函数值f(xi)
% f1a为f在左端点的一阶导数值
% f1b为f在右端点的一阶导数值
syms X;
n = length(x);
h = diff(x);
ff = diff(f)./h; %ff表示三弯矩方法中的f(x(i),x(i+1));
mu(1,1:n-2) = h(1,1:n-2)./(h(1,1:n-2)+h(1,2:n-1));
lambda(1,2:n-1) = 1-mu(1,1:n-2);
mu(1,n-1) = 1; %μ最后一位为1
lambda(1,1) = 1; %λ第一位为1
d = 6./(h(1,1:n-2)+h(1,2:n-1)).*(ff(1,2:n-1)-ff(1,1:n-2));
d0 = 6/h(1)*(ff(1)-f1a); %计算d0
dn = 6/h(n-1)*(f1b-ff(n-1)); %计算dn
d = [d0,d,dn];
r(1) = 2;
z(1) = d(1);
for i=2:n %追
l(i) = mu(i-1)/r(i-1);
r(i) = 2-l(i)*lambda(i-1);
z(i) = d(i)-l(i)*z(i-1);
end
m(n)=z(n)/r(n);
for j=n-1:-1:1 %赶
m(j) = (z(j)-lambda(j)*m(j+1))/r(i);
end
for j=2:1:n %讲求得的m带入书80页式4.7.2,求得S(x)
S(j-1)=vpa(((x(j)-X).^3)*m(j-1)/6/h(j-1)+((X-x(j-1)).^3)*m(j)/6/h(j-1)+(f(j-1)-m(j-1)*(h(j-1).^2)/6)*(x(j)-X)/h(j-1)+(f(j)-m(j)*(h(j-1).^2)/6)*(X-x(j-1))/h(j-1),2);
fprintf('s(%d)=%1.2f',j-1);
pretty(collect(S(j-1))) %输出方程形式
end
end