I’m working with a simple RBC model specified specified in levels so that Dynare linearizes it.
I already know the solutions for capital (K) and labour (L).
The output (Y) equation is not specified in the set of model equations as Y was substituted out.

How can I obtain the impulse response of Y without adding it to the model equations since I don’t need the first order approximation to Y as I already know K and L?

Y is computed by using the production function Y= K^(1-alpha)*L^(alpha).

From your model you get the irf for K and L and the corresponding steady state values save in oo_.steady_state().
Then just use these values to calculate the irf for Y:

Y_irf = (K_eps_YOURSHOCK + oo_.steady_state(HERE THE NUMBER WHERE DYNARE HAS SAVED THE SS FOR K).ones(size(K_eps_YOURSHOCK)))^(1-alpha) (L_eps_YOURSHOCK + oo_.steady_state(HERE THE NUMBER WHERE DYNARE HAS SAVED THE SS FOR L).*ones(size(L_eps_YOURSHOCK)))^(alpha)

Seems there is an error in the expression: should coefficients be .^(1-alpha) and .^(alpha)?

Is the product of the two terms involving K and L conformable to give Y_irf of dimension size(K_eps_YOURSHOCK) x 1?

Y_irf = (K_eps_YOURSHOCK + oo_.steady_state(HERE THE NUMBER WHERE DYNARE HAS SAVED THE SS FOR K).ones(size(K_eps_YOURSHOCK)))^(1-alpha) (L_eps_YOURSHOCK + oo_.steady_state(HERE THE NUMBER WHERE DYNARE HAS SAVED THE SS FOR L).*ones(size(L_eps_YOURSHOCK)))^(alpha)

If Y=AK^(1-alpha)*L^(alpha), where A is the AR(1) shock, then one need to multiply the expression above with (A_eps_YOURSHOCK +oo_.steady_state(HERE THE NUMBER WHERE DYNARE HAS SAVED THE SS FOR A))?