let {i in A} a[i] := 0;
let a[1] := -1000;
let a[2] :=  2000;
let a[3] := -1000;

let {i in N} c[i] := F[i] - K[i];

option filter_options "maxiter=20000";
option ipopt_options 'iprint=0 imaxiter=20000 dtol=1e-6';
option snopt_options "iterations_limit=500000 Major_iterations= 50000";



solve;

display a;
display {i in N} (K[i], c[i], util[i], sum {j in A} a[j]*K[i]^(j-1), sum {j in A} a[j]*(F[i]-c[i])^(j-1), Vgrad.body[i], Vhess.body[i], BellmanEQ.body[i]);

display _solve_time ;  


# Check Bellman Equation errors on uniform grid
# L-infinity norm
display max {i in UN} abs(sum {j in A} a[j]*(U[i]^(j-1)) - ((sum {k in A: k >= 2} a[k]*(k-1)*(U[i]^(k-2))/(1+(1-beta)/beta*U[i]^(alpha-1)))^(1/gamma))^(gamma+1)/(gamma+1) - beta*(sum {j in A} a[j]*(UF[i]-(sum {k in A: k >= 2} a[k]*(k-1)*(U[i]^(k-2))/(1+(1-beta)/beta*U[i]^(alpha-1)))^(1/gamma))^(j-1)));

# L-2 norm
display (sum {i in UN} (sum {j in A} a[j]*(U[i]^(j-1)) - ((sum {k in A: k >= 2} a[k]*(k-1)*(U[i]^(k-2))/(1+(1-beta)/beta*U[i]^(alpha-1)))^(1/gamma))^(gamma+1)/(gamma+1) - beta*(sum {j in A} a[j]*(UF[i]-(sum {k in A: k >= 2} a[k]*(k-1)*(U[i]^(k-2))/(1+(1-beta)/beta*U[i]^(alpha-1)))^(1/gamma))^(j-1)))^2)/uk;