E.g. y'(x) + ay(x) = f(x) where a is a constant.

 

I know we basically get exp(ax), then equation becomes [exp(ax)y(x)]' = exp(ax)f(x). But since f(x) is abstract/not defined, I don't know how I would not get stuck with a recursive integral on RHS after trying to integrate.

Anyone know how to evaluate integral of exp(ax)f(x)?

You have a nonhomogeneous linear differential equation.

The standard routines can be found here

https://www.google.co.uk/search?q=nonhomogeneous+first+order+linear+differential+equation&ie=utf-8&oe=utf-8&client=firefox-b&gfe_rd=cr&dcr=0&ei=yT8wWvzGJ6zP8Ae7g4vAAQ

Well, I see the procedure is different. However, looking up the particular solution to use would be easier if I actually knew what f(x) was. I am basically trying to solve for neutron population with a source term but the source function isn't known. I would hazard a guess as to exponential nature which seems most accurate but is there any reason to believe it would not be a polynomial?

Since f(x) is unknown, have you thought about examining the direction field, and the orthogonal trajectories to it, to find out what form it must have?

  Using "variation of parameters":  Given that y'+ ay= f(x), the "associated homogeneous equation" is y'+ ay= 0.  y'= dy/dx= -ay so dy/y= -a dx.  Integrating ln(y)= -ax+ C, y= C'e^-ax

Now we look for a solution to the entire equation of the form y(x)= u(x)e^-ax

Then y'= u'e^-ax- aue^-ax.  Putting that into the equation, u'e^-ax- aue^-ax+ aue^-ax= u'e^-ax= f(x).  u'(x)= e^axf(x) so u(x)= integral e^ax f(x) dx.  The general solution to the given differential equation is y(x)= C'e^-ax+ e^-ax(integral e^ax f(x)dx).