function dp2drho2 = eval_ppp(x,v,model, Te0, ne0, Ti0 , ni0)
%#codegen

% Tep, Tepp
tehat = x(model.te.xind);
tep0 = model.te.Lamp(1,:)*tehat;
tepp0 = model.te.Lampp(1,:)*tehat;

% ti is scaled version of ne.
% scaling of x coefficients is stored in v
scal = model.ti.Lam(1,:)*v(model.ti.vind,:);
tipp0 =  scal.*tepp0;
tip0 = scal.*tep0;
nehat = v(model.ne.vind);
nep0 = model.ne.Lamp(1,:)*nehat;
nepp0 = model.ne.Lampp(1,:)*nehat;
% ni is directly specified by v
nip0 = model.ni.Lamp(1,:)*v(model.ni.vind,:);
nipp0 = model.ni.Lampp(1,:)*v(model.ni.vind,:);


dp2drho2 = 1.602e-19*(tepp0*ne0+nepp0*Te0 + 2*tep0*nep0...
    + tipp0*ni0 + nipp0*Ti0 + 2*tip0*nip0);
    
return