Adding symbolic functions and manipulating results of DSolve

Question

I want to add a symbolic function and its second derivative. The symbolic function is given by DSolve:

sol = DSolve[y''[x] + w^2 y[x] == A Cos[q x], y[x], x]
(* {{y[x] -> -((A Cos[q x])/(q^2 - w^2)) + C[1] Cos[w x] + 
    C[2] Sin[w x]}} *)

p[x_, q_, w_, A_] = y[x] /. sol[[1]]
(* -((A Cos[q x])/(q^2 - w^2)) + C[1] Cos[w x] + C[2] Sin[w x] *)

g[x_, q_, w_, A_] = p''[x, q, w, A] + w^2 p[x, q, w, A]
(* {((y'')[x] -> (A q^2 Cos[q x])/(q^2 - w^2) - 
     w^2 C[1] Cos[w x] - w^2 C[2] Sin[w x]) + 
  w^2 (-((A Cos[q x])/(q^2 - w^2)) + C[1] Cos[w x] + C[2] Sin[w x])} *)

I want Mathematica to combine like terms in g[x,q,w,A] when displayed, and I want to be able to evaluate g[1,2,3,4] to get a single numeric value, like so:

N[p[1, 2, 3, 4]]
(* -0.332917 - 0.989992 C[1] + 0.14112 C[2] *)

Which is very different from the form of the result given by g[1,2,3,4], as it contains (y''[x])[1]-> ...:

g[1, 2, 3, 4]
(* {((y''[x])[1] -> -((16 Cos[2])/5) - 9 C[1] Cos[3] - 9 C[2] Sin[3]) + 9 ((4 Cos[2])/5 + C[1] Cos[3] + C[2] Sin[3])} *)

How do I make sure g works like a regular symbolic function?

Answer 1

$Version

(* "13.3.1 for Mac OS X ARM (64-bit) (July 24, 2023)" *)

Clear["Global`*"]

eqn = y''[x] + w^2 y[x] == A Cos[q x];

Ask for the solution as a pure function

sol = DSolve[eqn, y, x]

(* {{y -> Function[{x}, -((A*Cos[q*x])/(q^2 - w^2)) + 
           C[1]*Cos[w*x] + C[2]*Sin[w*x]]}} *)

Derivatives evaluate automatically,

eqn /. sol[[1]] // Simplify

(* True *)

p[x_, q_, w_, A_] = y[x] /. sol[[1]]

(* -((A*Cos[q*x])/(q^2 - w^2)) + 
   C[1]*Cos[w*x] + C[2]*Sin[w*x] *)

In the definition of g, p''[x, q, w, A] is undefined since you have not specified which variable to use in the derivative. For the second derivative with repect to x, use Derivative[2, 0, 0, 0][p][x, q, w, A]

g[x_, q_, w_, A_] = 
 Derivative[2, 0, 0, 0][p][x, q, w, A] + w^2 p[x, q, w, A] // Simplify

(* A Cos[q x] *)

N[p[1, 2, 3, 4]]

(* -0.3329174692377139 - 0.9899924966004454*C[1] + 
   0.1411200080598672*C[2] *)

g[1, 2, 3, 4]

(* 4 Cos[2] *)