%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
% function quad(xx,h,para,'u')
%
% makoto nakajima, October 11,1999
%
% This routine culculate quadratic approximation of the function 'u'
% around steady state xs.
%
% This program is translated from GAUSS program of Gary Hansen
%
%
% locally use : i,j,n,q,z,hz
% inputs : u (function to be approximated)
% xx (point at which the function is approximated)
% h (a small positive number)
% para(parameters in f)
% output : r_matrix
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

function[r_matrix]=vquad(xx,h,para,u) 

n=length(xx);
q=zeros(n+1,n+1);
hz=abs(h*xx); 

for i=1:n
z(i,i)=hz(i);
if hz(i)==0
z(i,i)=h;
end
end 

i=1; 

while i<=n
q(i+1,i+1)=((feval(u,xx+z(:,i)',para)+feval(u,xx-z(:,i)',para))/2-feval(u,xx,para))/z(i,i)^2;
j=1;
while j<i;
q(i+1,j+1)=(feval(u,xx+z(:,i)'+z(:,j)',para)-feval(u,xx+z(:,i)'-z(:,j)',para)-feval(u,xx-z(:,i)'+z(:,j)',para)+feval(u,xx-z(:,i)'-z(:,j)',para))/(8*z(i,i)*z(j,j));
q(j+1,i+1)=q(i+1,j+1);
j=j+1;
end
i=i+1;
end 

i=1;
while i<=n
q(i+1,1)=(feval(u,xx+z(:,i)',para)-feval(u,xx-z(:,i)',para))/(2*z(i,i));
q(i+1,1)=(q(i+1,1)-2*(q(i+1,2:n+1)*xx'))/2;
q(1,i+1)=q(i+1,1);
i=i+1;
end
r_matrix=q;