proc quadp(xs,h,&u);
@Calculate Quadratic Approximation

     xs=vector containing steady state values.
     h=small positive constant
     u=proc containing function to be approximated @

local i,j,n,q,z,hz,u:proc;
n=rows(xs);
q=zeros(n+1,n+1);
z=zeros(n,n);
hz=h*ones(n,1);
z=diagrv(z,hz);
i=1;
do until i>n;
     q[i+1,i+1]=((u(xs+z[.,i])+u(xs-z[.,i]))/2-u(xs))/z[i,i]^2;
     j=1;
     do until j==i;
          q[i+1,j+1]=(u(xs+z[.,i]+z[.,j])-u(xs+z[.,i]-z[.,j])
                     -u(xs-z[.,i]+z[.,j])+u(xs-z[.,i]-z[.,j]))/
                     (8*z[i,i]*z[j,j]);
          q[j+1,i+1]=q[i+1,j+1];
          j=j+1;
     endo;
     i=i+1;
endo;
i=1;
do until i>n;
     q[i+1,1]=(u(xs+z[.,i])-u(xs-z[.,i]))/(2*z[i,i]);
     q[i+1,1]=(q[i+1,1]-2*(q[i+1,2:n+1]*xs))/2;
     q[1,i+1]=q[i+1,1];
     i=i+1;
endo;

retp(q);
endp;