%input of system of ordinary differential equations (as a vector)
 
funcODE=[ output_prm(k,t)-c-p*k-ddelta*k;...       %ODE: k
          c*((1/ssigma)*(mpk_prm(k,t)-rrho-ddelta))-p*m*k; %ODE: c
          m*(m-(mpk_prm(k,t)-ddelta+(1/ssigma)*(rrho+ddelta-mpk_prm(k,t))))];
      %z*((1/gamma-1)*(alpha*A_neutral*alpha^(-alpha)*(1-alpha)^(alpha-1)*k^(alpha-1)*ell^(1-alpha)-delta)+z-rho/gamma)   
      %z-(alpha*A_neutral*alpha^(-alpha)*(1-alpha)^(alpha-1)*k^(alpha-1)*ell^(1-alpha)-delta)*(1-1/gamma)-(1/gamma)*rho
         %first state, then control