clc
clear all
close all
%%A*y[n]+B*y[n-1]+C*y[n-2]=Am1*x[n]+Am2*x[n-1]+Am3*x[n+1]
%x[n]=Am*a^n
syms t n k k1 k2 ;
A=1;  %y[n]
B=-0.6;%y[n-1]
C=0;%y[n-2]
a=0.4;%a^n
i1=10;%y[-1]
i2=0;%y[-2]
Am=1
m1=n;
Am1=1;%Am1*x[n]
m2=n-1;
Am2=0;%Am2*x[n-1]
m3=n+1;
Am3=0;%Am3*x[n+1]

xt1=Am*a^n;

y1='A+B/t+C/(t^2)=0';
f=subs(sym(y1));
f=char(f);
z=solve(f);

if C~=0
    if z(1)~=z(2)&imag(z(1))==0
    yn=k1*(z(1)^n)+k2*(z(2)^n);    
elseif z(1)~=z(2)&imag(z(1))~=0
    omega=angle(double(z(1)));
    r=abs(z(1));
    yn=r^n*(k1*cos(n*omega)+k2*sin(n*omega));
elseif z(1)==z(2)&imag(z(1))==0
    yn=z(1)^n*(k1+k2*n);
elseif z(1)==z(2)&imag(z(1))~=0
    ;
end
elseif C==0
    yn=k*(z^n);
end

y2=sym('A*c*a^n+B*c*a^n/a+C*c*a^n/a^2');
f=subs(y2);
f1=sym('f/xt1=1');
f=subs(f1);
f=char(f);
c=solve(f);
yf=sym('c*xt1');
yf=subs(yf);
yf=char(yf);

y=yn+yf;
yzi=yn;

f1='y=i1';
fi1='yzi=i1';
n=-1;
f1=subs(sym(f1));
f1=subs(f1);
f1=char(f1);
fi1=subs(sym(fi1));
fi1=subs(fi1);
fi1=char(fi1);

f2='y=i2';
fi2='yzi=i2';
n=-2;
f2=subs(sym(f2));
f2=subs(f2);
f2=char(f2);
fi2=subs(sym(fi2));
fi2=subs(fi2);
fi2=char(fi2);

if C~=0
s=solve(f1,f2,'k1','k2');
s1=solve(fi1,fi2,'k1','k2');

k1=s.k1;
k2=s.k2;
ki1=s1.k1;
ki2=s1.k2;
else 
    k=solve(f1);
    ki1=solve(fi1);
end
syms n;
yn=subs(sym(yn));
yn=char(yn);
if C~=0
k1=ki1;
k2=ki2;
else k=ki1;
end
yzi=subs(sym(yzi));
yzi=char(yzi);
yzii=yzi;
y='yf+yn';
yzs='y-yzi';
y=subs(sym(y));
y=char(y);
yzs=subs(sym(yzs));
yzs=char(yzs);

yy=inline(y,'n');
yyzi=inline(yzi,'n');
yyzs=inline(yzs,'n');
yyf=inline(yf,'n');
yyn=inline(yn,'n');
y=yy;
yzi=yyzi;
yzs=yyzs;
yf=yyf;
yn=yyn;

yzi=yzi(m1)*'u(m1)';
yzs=Am1*yzs(m1)*'u(m1)'+Am2*yzs(m2)*'u(m2)'+Am3*yzs(m3)*'u(m3)';
y=subs(sym(y))
yzi=subs(sym(yzi))
yzs=subs(sym(yzs))


yzs1=Am1*yyzs(m1)+Am3*yyzs(m3);
yzs2=yzs1+Am2*yyzs(m2);
n=0:0.1:0.9;
yzs=subs(sym(yzs1));
m=zeros(1,91);
ymzs1=[yzs,m];
n=1:0.1:10;
yzs=subs(sym(yzs2));
m=zeros(1,10);
ymzs2=[m,yzs];
yzs=ymzs1+ymzs2;
n=0:0.1:10;
yzi=subs(sym(yzii));
yf=subs(sym(yf));
figure
subplot(2,3,1)
plot(n,y);
title('全響應')
subplot(2,3,2)
plot(n,yf.*n./n);
title('受迫響應')
subplot(2,3,3)
plot(n,yn.*n./n);
title('自然響應')
subplot(2,3,4)
plot(n,yzi.*n./n);
title('零輸入響應')
subplot(2,3,5)
plot(n,yzs.*n./n);
title('零狀態(tài)響應')