PadeApprox_Ration.m

function [y,cx,H,coffsab,Rnx]=PadeApprox_Ration(f,n,m)
% pade逼近计算程序,
syms x

y=taylor(f,x,'Order',n+m+1);
coffsy0=sym2poly(y);p=length(coffsy0);coffsy=[zeros(1,n+m+1-p),coffsy0];
cx=flip(coffsy);

H2=zeros(n+1,m);
for i=2:m+1
     r=m+2-i;
     H2(i:n+1,r)=-cx(1:n+2-i)';
end
H1=eye(n+1);
H3=zeros(m,n+1);
H4=zeros(m,m);
for i=1:m
    H4(i,:)=-cx(n-m+1+i:n+i);
end
H=[H1,H2;H3,H4];
coffsab=H\cx';%coffsab=[a0,a1,***an,bm,bm-1,...b1]

j=0:n;Funa=x.^j;%定义分子
k=0:m;Funb=x.^k;%定义分母
digits(5);
coffsab2=flip(coffsab(n+2:n+m+1));
coffsab3=[1;coffsab2];
Rnx=vpa((Funa*coffsab(1:n+1))/(Funb*coffsab3));