function [nfx,nfxp,nfy,nfyp,nfypyp,nfypy,nfypxp,nfypx,nfyyp,nfyy,nfyxp,nfyx,nfxpyp,nfxpy,nfxpxp,nfxpx,nfxyp,nfxy,nfxxp,nfxx] = neoclassical_model
%This program computes numerical first and second derivatives of the  function f for the simple neoclassical growth model described in section 2.1 of ``Solving Dynamic General Equilibrium Models Using a Second-Order Approximation to the Policy Function,'' by Stephanie Schmitt-Grohe and Martin Uribe, (2001). Unlike the example in section 2.1, here y and x are defined as log(c) and [log(k); log(A)] respectively.   The function f defines  the DSGE model: 
%  E_t f(yp,y,xp,x) =0. 
%
%Inputs: none
%
%Output: Numerical first and second derivatives of f
%
%Calls: anal_deriv.m num_eval.m neoclassical_model_ss.m
%
%(c) Stephanie Schmitt-Grohe and Martin Uribe
%Date July 17, 2001

%Define parameters
syms 	SIG DELLTA ALFA BETTA RHO

%Define variables 
syms c cp k kp a ap 

%Write equations ei, i=1:3
e1 = exp(c) + exp(kp) - (1-DELLTA) * exp(k) - exp(a) * exp(k)^ALFA;
e2 = exp(c)^(-SIG) - BETTA * exp(cp)^(-SIG) * (exp(ap) * ALFA * exp(kp)^(ALFA-1) + 1 - DELLTA);
e3 = ap - RHO * a;

%Create function f
f = [e1;e2;e3];

% Define the vector of controls, y, and states, x
x = [k a];
y = c;
xp = [kp ap];
yp = cp;
nx = length(x);
ny = length(y);

%Compute analytical derivatives of f
[fx,fxp,fy,fyp,fypyp,fypy,fypxp,fypx,fyyp,fyy,fyxp,fyx,fxpyp,fxpy,fxpxp,fxpx,fxyp,fxy,fxxp,fxx]=anal_deriv(f,x,y,xp,yp);

%Numerical Evaluation
%Steady State and Parameter Values
[SIG,DELLTA,ALFA,BETTA,RHO,c,cp,k,kp,a,ap]=neoclassical_model_ss;
%Order of approximation desired 
approx = 2;
%Obtain numerical derivatives of f
num_eval