function f = foc(P,P0,Y,T,g,d)
for i=1:T-1
    if i==1
        f(i)=g*Y(i)+P0-(1+g+d)*P(i)+d*P(i+1);
    else
        f(i)=g*Y(i)+P(i-1)-(1+g+d)*P(i)+d*P(i+1);
    end
end
    f(T)=g*Y(T)+P(T-1)-(1+g+d)*P(T)+d*Y(T);