function [f,g] = rody(x);
%RODY     The objective of Jones-Manuelli-Rossi, JPE 1993 for the simple case.

T          = 50;
alpha      = .36;
beta       = .98;
delk       = .1;
delh       = .1;
tauk       = .21;
tauh       = .31;
gy         = .2;
ty         = .0726;
h0         = 1.0;
k0         = alpha*(1-tauk)*h0/((1-alpha)*(1-tauh));
k0         = alpha*h0/(1-alpha);
sigma      = 2.0;
A          = .43;
%gamma     = (beta*(A*((1-alpha)*(1-tauh))^(1-alpha)*(alpha*(1-tauk))^alpha ...
%            +1-delk))^(1/sigma);
gamma      = (beta*(A*(1-alpha)^(1-alpha)*alpha^alpha+1-delk))^(1/sigma);
hk         = h0/k0;
ck         = (1-gy)*A*hk^(1-alpha)-delk-delh*hk;
W0         = (1-delk+(1-tauk)*alpha*A*hk^(1-alpha))*k0 ...
            +(1-delh+(1-tauh)*(1-alpha)*A*hk^(-alpha))*h0;
k          = [k0;x(1:T+1)]; 
h          = [h0;x(T+2:2*T+2)];
gstar      = beta^(T+1)*(ck-ty*A*hk^(1-alpha))*ck^(-sigma)* ...
             k(T+2)^(1-sigma)/(1-beta*gamma^(1-sigma));
fstar      = beta^(T+1)*(ck*k(T+2))^(1-sigma)/(1-sigma)/ ...
             (1-beta*gamma^(1-sigma));

% u(c) = c^(1-sigma)/(1-sigma),  f(k,h) = A k^alpha*h^(1-alpha)

% 0,1,...T,  T+1 <-- in article
% 1,2,   T+1,T+2 <-- use

y              = A*k.^alpha.*h.^(1-alpha);
c              = (1-gy)*y(1:T+1)-k(2:T+2)+(1-delk)*k(1:T+1)-h(2:T+2)+ ...
                 (1-delh)*h(1:T+1);
%ink           = -(1-delk)*k(1:T)+k(2:T+1);
%inh           = -(1-delh)*h(1:T)+h(2:T+1);

W0             = W0*c(1)^(-sigma);

%[sum( (beta.^[0:T]').*(c(1:T+1)-ty*y(1:T+1)).*c(1:T+1).^(-sigma) ),gstar,W0]

g(1)           = sum( (beta.^[0:T]').*(c(1:T+1)-ty*y(1:T+1)).* ...
                 c(1:T+1).^(-sigma) ) + gstar - W0;
g(2)           = h(T+2)/k(T+2)-hk;
g(3:T+2)       = (1-delk)*k(1:T)-k(2:T+1);
g(T+3:2*T+2)   = (1-delh)*h(1:T)-h(2:T+1);
g(2*T+3:3*T+2) = c(1:T)-c(2:T+1);
g(3*T+3:4*T+3) = -c(1:T+1);
g              = g(:);
f              = -sum((beta.^[0:T]').*c.^(1-sigma)./(1-sigma)) - fstar;

plot(0:T+1,log(h),0:T+1,log(k),0:T,log(c)),pause(.1)

