% Precise Integration Method
clear;
D=[0.5,0;0,1];
B=[-6,2;2,-4];
f0=[0;0;0;10];
f1=zeros(size(f0));
H=[A,D;B,C];
I=eye(size(H));
iH=inv(H);
step=[2,0.5,0.1]; % different step size
Ta=H*dt+(H*dt)^2*(I+(H*dt)/3+(H*dt)^2/12)/2;
t(:,iter)=step(jj)*(iter-1);
vk=T*(vk+iH*(f0+iH*f1))-iH*(f0+iH*f1+f1*step(jj));
plot(t(1:tf/step(jj)),v,str(jj));
end