Hi Michael
Thank you for the detailed response.
I managed to develop the below code which fits my needs really well. I am just posting it here in case someone ever has a similar query.
from wolframclient.evaluation import WolframLanguageSession
from wolframclient.language import wl, wlexpr
import matplotlib.pyplot as plt
import numpy as np
session = WolframLanguageSession()
wolfram_code = """
s = NDSolve[{y'[x] == y[x] Cos[x + y[x]], y[0] == 1}, y, {x, 0, 30}];
Table[{x, y[x] /. s}, {x, 0, 30, 0.1}]
"""
# Evaluate the Wolfram code
result = session.evaluate(wlexpr(wolfram_code))
# Terminate the Wolfram session
session.terminate()
# Extract x and y values from the result
x_vals = [pair[0] for pair in result]
y_vals = [pair[1] for pair in result]
# Plot the solution using matplotlib
plt.plot(x_vals, y_vals, label='y(x)')
plt.xlabel('x')
plt.ylabel('y(x)')
plt.title('Solution of y\'[x] = y[x] Cos[x + y[x]]')
plt.grid(True)
plt.legend()
# Show the plot
plt.show()