n=zeros(101,1);
C1=zeros(101,1);
C2=zeros(101,1);
C3=zeros(101,1);
C4=zeros(101,1);
C5=zeros(101,1);
C6=zeros(101,1);
n(1)=1;
l=0.00002;
C1(1)=0.000266/l/0.0127;
C2(1)=0.001491/l/0.0317;
C3(1)=0.001316/l/0.115;
C4(1)=0.002849/l/0.311;
C5(1)=0.000896/l/1.4;
C6(1)=0.000182/l/3.87;
h=0.01;
rct=0.003;
F=[(rct-0.007)/l,0.0127,0.0317,0.115,0.311,1.4,3.87;
0.000266/l,-0.0127,0,0,0,0,0;
0.001491/l,0,-0.0317,0,0,0,0;
0.001316/l,0,0,-0.115,0,0,0;
0.002849/l,0,0,0,-0.311,0,0;
0.000896/l,0,0,0,0,-1.4,0;
0.000182/l,0,0,0,0,0,-3.87];
A=[1,1,1,1,1,1,1;
0,1,2,3,4,5,6;
0,0,1,3,6,10,15;
0,0,0,1,4,10,20;
0,0,0,0,1,5,15;
0,0,0,0,0,1,6;
0,0,0,0,0,0,1];
I=ones(7);
L=[720/1764;1;1624/1764;735/1764;175/1764;21/1764;1/1764];
for i=1:1:100
y=[n(i);C1(i);C2(i);C3(i);C4(i);C5(i);C6(i)];
y1=h.*(F*y);
y2=(h^2/factorial(2)).*(F^2*y);
y3=(h^3/factorial(3)).*(F^3*y);
y4=(h^4/factorial(4)).*(F^4*y);
y5=(h^5/factorial(5)).*(F^5*y);
y6=(h^6/factorial(6)).*(F^6*y);
Z=[y,y1,y2,y3,y4,y5,y6]';
Z0=A*Z;
G=-Z0(2,:)'+h.*F*Z0(1,:)';
N=(L(2).*I-(h*L(1)).*F')^(-1);
Z1=Z0+dot(L,G').*N;
G=-Z1(2,:)'+h.*F*Z1(1,:)';
N=(L(2).*I-(h*L(1)).*F')^(-1);
Z2=Z1+dot(L,G').*N;
G=-Z2(2,:)'+h.*F*Z2(1,:)';
N=(L(2).*I-(h*L(1)).*F')^(-1);
Z3=Z2+dot(L,G').*N;
n(i+1)=Z0(1,1);
C1(i+1)=Z0(1,2);
C2(i+1)=Z0(1,3);
C3(i+1)=Z0(1,4);
C4(i+1)=Z0(1,5);
C5(i+1)=Z0(1,6);
C6(i+1)=Z0(1,7);
end
t=0:0.01:1;
plot(t,n)
