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);